{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# 纯规范场（淬火近似）杂化 Monte Carlo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 0 Markov chain Monte Carlo (MCMC)\n",
    "\n",
    "根据格点上的路径积分形式，首先完成费米子部分积分\n",
    "$$\n",
    "\\mathcal{Z}=\\int\\mathcal{D}U_\\mu e^{-S_g(U_\\mu)}\\det\\left[\\mathcal{M}(U_\\mu)\\right]\\\\\n",
    "\\langle\\mathcal{O}\\rangle _\\text{fermion}=\\det\\left[\\mathcal{M}(U_\\mu)\\right]^{-1}\\int\\mathcal{D}\\bar\\psi\\mathcal{D}\\psi e^{-\\bar\\psi\\mathcal{M}(U_\\mu)\\psi}\\mathcal{O}\\\\\n",
    "\\langle\\mathcal{O}\\rangle=\\frac{1}{\\mathcal{Z}}\\int\\mathcal{D}U_\\mu e^{-S_g(U_\\mu)}\\det\\left[\\mathcal{M}(U_\\mu)\\right]\\langle\\mathcal{O}\\rangle _\\text{fermion}\n",
    "$$\n",
    "其中 $\\langle\\mathcal{O}\\rangle _\\text{fermion}$ 部分可以通过在某个规范场组态（configuration）上进行测量得到。\n",
    "\n",
    "现在考虑规范场积分步骤。显然 $U_\\mu$ 的维数过于巨大，我们只能采用 Monte Carlo 方法进行数值积分。重点抽样可以极大程度增加积分的精度，因此我们希望有一些规范场满足分布\n",
    "$$\n",
    "\\left\\{U_\\mu\\right\\}\\sim\\frac{1}{\\mathcal{Z}}e^{-S_g(U_\\mu)}\\det\\left[\\mathcal{M}(U_\\mu)\\right]\n",
    "$$\n",
    "这样，我们只需要进行简单的对组态平均即可完成积分\n",
    "$$\n",
    "\\frac{1}{\\mathcal{Z}}\\int\\mathcal{D}U_\\mu e^{-S_g(U_\\mu)}\\det\\left[\\mathcal{M}(U_\\mu)\\right]\\langle\\mathcal{O}\\rangle _\\text{fermion}\\approx\\frac{1}{N}\\sum_{\\left\\{U_\\mu\\right\\}}\\langle\\mathcal{O}\\rangle _\\text{fermion}\n",
    "$$\n",
    "\n",
    "对于纯规范场（淬火近似）的组态来说，有 $\\det\\left[\\mathcal{M}(U_\\mu)\\right]=1$，上述分布退化到\n",
    "$$\n",
    "\\mathcal{Z}=\\int\\mathcal{D}U_\\mu e^{-S_g(U_\\mu)},\\;\\left\\{U_\\mu\\right\\}\\sim\\frac{1}{\\mathcal{Z}}e^{-S_g(U_\\mu)}\n",
    "$$\n",
    "现在分布形式与统计力学中 Boltzmann 分布十分类似。这样我们可以通过 Markov 链 Monte Carlo 方法获得一些符合目标分布的规范场组态。\n",
    "\n",
    "需要注意的是，通过 MCMC 径迹（trajectory）采样获得的样本之间存在自相关（autocorrelation），直接进行样本平均将会低估误差，需要额外的步骤降低影响。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 1 初始化"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 1.1 导入需要的 Python 包\n",
    "\n",
    "`exp` 和 `random` 用于判断是否接受新的规范场\n",
    "\n",
    "`perf_counter` 函数用于记时\n",
    "\n",
    "`pyquda` 包中的函数用于使用 GPU 加速格点 QCD 计算\n",
    "\n",
    "`hmc` 包来自 `hmc.py`，主要是对 PyQUDA 函数的封装"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "from math import exp\n",
    "from random import random\n",
    "from time import perf_counter\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from pyquda_utils import core, io\n",
    "\n",
    "from hmc import HMC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 1.2 初始化 PyQUDA 并设置格子大小\n",
    "\n",
    "调用 `init()` 初始化 PyQUDA，`resource_path` 指向存储调优参数的文件夹。\n",
    "\n",
    "获取一个 `LatticeInfo` 实例用于之后的计算，这里我们设置了一个 $16^3\\times32$ 的格子。\n",
    "\n",
    "获取一个 `HMC` 实例并初始化。HMC 将 QUDA 函数包装为更易用的形式。\n",
    "\n",
    "获取一个 `LatticeGauge` 实例用于存取每一步的规范场。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "PyQUDA INFO: Using the grid size [1, 1, 1, 1]\n",
      "PyQUDA INFO: Using QUDA_RESOURCE_PATH=.cache\n",
      "PyQUDA INFO: Using CUDA backend cupy\n",
      "Disabling GPU-Direct RDMA access\n",
      "Enabling peer-to-peer copy engine and direct load/store access\n",
      "QUDA 1.1.0 (git 1.1.0-341431726-sm_60)\n",
      "CUDA Driver version = 12020\n",
      "CUDA Runtime version = 11080\n",
      "Graphic driver version = 535.183.01\n",
      "Found device 0: Tesla P100-PCIE-16GB\n",
      "Found device 1: Tesla P100-PCIE-16GB\n",
      "Found device 2: Tesla P100-PCIE-16GB\n",
      "Found device 3: Tesla P100-PCIE-16GB\n",
      "Using device 0: Tesla P100-PCIE-16GB\n",
      "Initializing monitoring on device 0: Tesla P100-PCIE-16GB\n",
      "WARNING: Data reordering done on GPU (set with QUDA_REORDER_LOCATION=GPU/CPU)\n",
      "WARNING: The path \".cache\" specified by QUDA_RESOURCE_PATH does not exist or is not a directory.\n",
      "WARNING: Caching of tuned parameters will be disabled\n",
      "WARNING: Cache file not found.  All kernels will be re-tuned (if tuning is enabled).\n",
      "WARNING: Using device memory pool allocator\n",
      "WARNING: Using pinned memory pool allocator\n",
      "cublasCreated successfully\n"
     ]
    }
   ],
   "source": [
    "core.init(resource_path=\".cache\")\n",
    "latt_info = core.LatticeInfo([16, 16, 16, 32])\n",
    "\n",
    "hmc = HMC(latt_info)\n",
    "hmc.initialize()\n",
    "\n",
    "gauge = core.LatticeGauge(latt_info)\n",
    "Nc = latt_info.Nc"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 1.3 设置规范作用量参数"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 1.3.0 分子动力学\n",
    "\n",
    "首先我们有规范场 $U_\\mu$，引入共轭动量 $\\pi_\\mu$，考虑 Hamiltonian\n",
    "$$\n",
    "\\mathcal{H}=\\frac{1}{2}\\sum_{x,\\mu}\\mathrm{Tr}\\left[\\pi_\\mu(x)\\pi_\\mu(x)\\right]+\\frac{\\beta}{N_c}\\sum_{x,\\mu,\\nu\\gt\\mu}\\mathfrak{Re}\\mathrm{Tr}\\left[\\mathbb{I}_{N_c}-U_\\mu(x)A_{\\mu,\\nu}(x)\\right]\n",
    "$$\n",
    "共轭动量 $\\pi_\\mu(x)\\in\\mathfrak{su}(3)$，为无迹 Hermitian 矩阵。\n",
    "\n",
    "有 Hamiltonian 正则方程\n",
    "$$\n",
    "\\frac{\\mathrm{d}}{\\mathrm{d}t}\\pi_\\mu(x)=-\\frac{\\partial S_g(U_\\mu)}{\\partial U_\\mu(x)},\\;\\frac{\\mathrm{d}}{\\mathrm{d}t}U_\\mu(x)=i\\pi_\\mu(x)U_\\mu(x)\n",
    "$$\n",
    "\n",
    "我们希望通过分子动力学过程使得\n",
    "$$\n",
    "\\{\\pi_\\mu,U_\\mu\\}\\sim e^{-\\mathcal{H}}=e^{-\\frac{1}{2}\\sum_{x,\\mu}\\mathrm{Tr}\\left[\\pi_\\mu(x)\\pi_\\mu(x)\\right]}e^{-\\frac{\\beta}{N_c}\\sum_{x,\\mu,\\nu\\gt\\mu}\\mathfrak{Re}\\mathrm{Tr}\\left[\\mathbb{I}_{N_c}-U_\\mu(x)A_{\\mu,\\nu}(x)\\right]}\n",
    "$$\n",
    "显然 $\\pi_\\mu$ 和 $U_\\mu$ 的分布是互相独立的。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 1.3.1 Wilson 规范作用量\n",
    "\n",
    "$$\n",
    "S_g(U_\\mu)=\\frac{\\beta}{N_c}\\sum_{x,\\mu,\\nu\\gt\\mu}\\mathfrak{Re}\\mathrm{Tr}\\left[\\mathbb{I}_{N_c}-U_\\mu(x)U_\\nu(x+\\hat{\\mu})U^\\dagger_\\mu(x+\\hat{\\nu})U^\\dagger_\\nu(x)\\right]\n",
    "$$\n",
    "\n",
    "![plaquette](plaquette.drawio.svg)\n",
    "\n",
    "对于一个由规范链接收尾相连组成的圈，总可以使用方向矢量来表示。例如上图的圈可以写成 $\\mathrm{Tr}\\left[U_\\mu(x)U_\\nu(x+\\hat{\\mu})U^\\dagger_\\mu(x+\\hat{\\nu})U^\\dagger_\\nu(x)\\right]\\sim\\left[\\mu,\\nu,-\\mu,-\\nu\\right]\\text{ at }x$。\n",
    "\n",
    "忽略常数项，Wilson 规范作用量记为\n",
    "$$\n",
    "S_g(U_\\mu)=-\\frac{\\beta}{N_c}\\sum_{x,\\mu,\\nu\\gt\\mu}\\mathfrak{Re}\\left[\\mu,\\nu,-\\mu,-\\nu\\right]\n",
    "$$\n",
    "\n",
    "展开 $\\sum_{\\nu\\gt\\mu}$\n",
    "$$\n",
    "\\begin{aligned}\n",
    "S_g(U_\\mu)=-\\frac{\\beta}{N_c}\\sum_{x}\\mathfrak{Re}\\left\\{\\left[x,y,-x,-y\\right]+\\left[x,z,-x,-z\\right]+\\left[y,z,-y,-z\\right]\\right.\\\\\n",
    "\\left.+\\left[x,t,-x,-t\\right]+\\left[y,t,-y,-t\\right]+\\left[z,t,-z,-t\\right]\\right\\}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "我们使用 $0,1,2,3$ 来标记 $x,y,z,t$ 四个方向，$4,5,6,7$ 来标记相应的反方向\n",
    "$$\n",
    "\\begin{aligned}\n",
    "S_g(U_\\mu)=-\\frac{\\beta}{N_c}\\sum_{x}\\mathfrak{Re}\\left\\{\\left[0,1,4,5\\right]+\\left[0,2,4,6\\right]+\\left[1,2,5,6\\right]\\right.\\\\\n",
    "\\left.+\\left[0,3,4,7\\right]+\\left[1,3,5,7\\right]+\\left[2,3,6,7\\right]\\right\\}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "PyQUDA 会处理取实部操作，对于特定坐标 $x$，最终输入的作用量为\n",
    "$$\n",
    "\\begin{aligned}\n",
    "S'_g(U_\\mu)=-\\frac{\\beta}{N_c}\\left\\{\\left[0,1,4,5\\right]+\\left[0,2,4,6\\right]+\\left[1,2,5,6\\right]\\right.\\\\\n",
    "\\left.+\\left[0,3,4,7\\right]+\\left[1,3,5,7\\right]+\\left[2,3,6,7\\right]\\right\\}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\beta=\\frac{2N_c}{g^2}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "u_0 = 0.855453\n",
    "beta = 6.20\n",
    "input_path = [\n",
    "    [0, 1, 4, 5],\n",
    "    [0, 2, 4, 6],\n",
    "    [1, 2, 5, 6],\n",
    "    [0, 3, 4, 7],\n",
    "    [1, 3, 5, 7],\n",
    "    [2, 3, 6, 7],\n",
    "]\n",
    "input_coeff = [\n",
    "    -beta / Nc,\n",
    "    -beta / Nc,\n",
    "    -beta / Nc,\n",
    "    -beta / Nc,\n",
    "    -beta / Nc,\n",
    "    -beta / Nc,\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### 1.3.2 Wilson 规范力\n",
    "\n",
    "根据 Hamiltonian 正则方程\n",
    "$$\n",
    "\\delta U_\\mu(x)=i\\pi_\\mu(x)U_\\mu(x)\\delta t\n",
    "$$\n",
    "\n",
    "分子动力学积分需要 Hamiltonian（总能量）保持不变\n",
    "$$\n",
    "\\delta H=0\n",
    "$$\n",
    "\n",
    "那么对某个特定位置 $x$ 的部分变分\n",
    "$$\n",
    "\\begin{aligned}\n",
    "0=&\\sum_{\\mu}\\mathrm{Tr}\\left[\\delta\\pi_\\mu(x)\\pi_\\mu(x)\\right]\\\\\n",
    "&-\\frac{\\beta}{N_c}\\sum_{\\mu,\\nu\\ne\\mu}\\mathfrak{Re}\\mathrm{Tr}\\left[\\delta U_\\mu(x)A_{\\mu,\\nu}(x)+\\delta U_\\mu(x)A_{\\mu,-\\nu}(x)\\right]\\delta t\\\\\n",
    "=&\\sum_{x,\\mu}\\mathrm{Tr}\\left[\\delta\\pi_\\mu(x)\\pi_\\mu(x)\\right]\\\\\n",
    "&-\\frac{\\beta}{N_c}\\sum_{x,\\mu,\\nu\\ne\\mu}\\mathfrak{Re}\\mathrm{Tr}\\left[i\\pi_\\mu(x)U_\\mu(x)A_{\\mu,\\nu}(x)+i\\pi_\\mu(x)U_\\mu(x)A_{\\mu,-\\nu}(x)\\right]\\delta t\\\\\n",
    "=&\\sum_{x,\\mu}\\mathrm{Tr}\\left[\\delta\\pi_\\mu(x)\\pi_\\mu(x)\\right]\\\\\n",
    "&-\\frac{\\beta}{N_c}\\sum_{x,\\mu,\\nu\\ne\\mu}\\mathrm{Tr}\\left\\{i\\pi_\\mu(x)\\left[U_\\mu(x)A_{\\mu,\\nu}(x)+U_\\mu(x)A_{\\mu,-\\nu}(x)\\right]_\\text{Traceless,Anti-Hermitian}\\right\\}\\delta t\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "&\\mathfrak{Re}\\mathrm{Tr}\\left[i\\pi_\\mu(x)U_\\mu(x)A_{\\mu,\\nu}(x)\\right]\\\\\n",
    "=&\\frac{1}{2}\\mathrm{Tr}\\left[i\\pi_\\mu(x)U_\\mu(x)A_{\\mu,\\nu}(x)-iA^\\dagger_{\\mu,\\nu}U^\\dagger_\\mu(x)\\pi_\\mu(x)\\right]\\\\\n",
    "=&\\mathrm{Tr}\\left\\{i\\pi_\\mu(x)\\left[U_\\mu(x)A_{\\mu,\\nu}(x)\\right]_\\text{Anti-Hermitian}\\right\\}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "![plaquette](plaquette.drawio.svg)\n",
    "\n",
    "上式应当对任意的 $\\pi_\\mu(x)$ 成立\n",
    "$$\n",
    "\\delta\\pi_\\mu(x)=\\frac{\\beta}{N_c}\\sum_{\\nu\\ne\\mu}i\\left[U_\\mu(x)A_{\\mu,\\nu}(x)+U_\\mu(x)A_{\\mu,-\\nu}(x)\\right]_\\text{Traceless,Anti-Hermitian}\\delta t\n",
    "$$\n",
    "\n",
    "QUDA 中动量场实际为 $i\\pi_\\mu$，对规范力添加了一个负号并应用无迹反 Hermitian 操作\n",
    "$$\n",
    "F'_g(x,U_\\mu)=\\frac{\\beta}{N_c}\\sum_{\\nu\\ne\\mu}\\left\\{[\\mu,\\nu,-\\mu,-\\nu]+[\\mu,-\\nu,-\\mu,\\nu]\\right\\}\n",
    "$$\n",
    "\n",
    "以 $\\mu=x$ 为例，按照之前约定使用数字代表方向，对于特定坐标 $x$，最终输入的规范力为\n",
    "$$\n",
    "\\begin{aligned}\n",
    "F'_g(U_\\mu)=\\frac{\\beta}{N_c}\\left\\{[0,1,4,5]+[0,5,4,1]+[0,2,4,6]\\right.\\\\\n",
    "\\left.+[0,6,4,2]+[0,3,4,7]+[0,7,4,3]\\right\\}\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "input_path2 = [\n",
    "    [\n",
    "        [0, 1, 4, 5],\n",
    "        [0, 5, 4, 1],\n",
    "        [0, 2, 4, 6],\n",
    "        [0, 6, 4, 2],\n",
    "        [0, 3, 4, 7],\n",
    "        [0, 7, 4, 3],\n",
    "    ],\n",
    "    [\n",
    "        [1, 4, 5, 0],\n",
    "        [1, 0, 5, 4],\n",
    "        [1, 2, 5, 6],\n",
    "        [1, 6, 5, 2],\n",
    "        [1, 3, 5, 7],\n",
    "        [1, 7, 5, 3],\n",
    "    ],\n",
    "    [\n",
    "        [2, 4, 6, 0],\n",
    "        [2, 0, 6, 4],\n",
    "        [2, 5, 6, 1],\n",
    "        [2, 1, 6, 5],\n",
    "        [2, 3, 6, 7],\n",
    "        [2, 7, 6, 3],\n",
    "    ],\n",
    "    [\n",
    "        [3, 4, 7, 0],\n",
    "        [3, 0, 7, 4],\n",
    "        [3, 5, 7, 1],\n",
    "        [3, 1, 7, 5],\n",
    "        [3, 6, 7, 2],\n",
    "        [3, 2, 7, 6],\n",
    "    ],\n",
    "]\n",
    "input_coeff2 = [\n",
    "    beta / Nc,\n",
    "    beta / Nc,\n",
    "    beta / Nc,\n",
    "    beta / Nc,\n",
    "    beta / Nc,\n",
    "    beta / Nc,\n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 1.4 设置 Markov 链（径迹）参数\n",
    "\n",
    "Markov 链从 0 开始（单位场），总共运行 2000 步，前 500 步不考虑 Metropolis 算法而是直接接受，使得系统达到热平衡。\n",
    "\n",
    "每一步的演化时间为 1.0，每次演化（分子动力学积分）分成 20 小步以降低积分误差。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "start = 0\n",
    "stop = 2000\n",
    "warm = 500\n",
    "save = 50\n",
    "t = 1.0\n",
    "n_steps = 20"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## 2 HMC"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 2.1 输出初始状态统计信息\n",
    "\n",
    "小方格（plaquette）作为最简单的可观测量，表征规范场势能。\n",
    "$$\n",
    "\\text{Plaquette}(U_\\mu)=\\frac{1}{6VN_c}\\sum_{x,\\mu,\\nu\\gt\\mu}\\mathfrak{Re}\\mathrm{Tr}\\left[U_\\mu(x)U_\\nu(x+\\hat{\\mu})U^\\dagger_\\mu(x+\\hat{\\nu})U^\\dagger_\\nu(x)\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Trajectory 0:\n",
      "plaquette = 1.0\n",
      "\n"
     ]
    }
   ],
   "source": [
    "plaquette = [hmc.plaquette()]\n",
    "print(\"\\n\" f\"Trajectory {start}:\\n\" f\"plaquette = {plaquette[-1]}\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 2.2 第 1 步开始\n",
    "\n",
    "计数器设置为 `start`。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trajectory 1:\n",
      "\n"
     ]
    }
   ],
   "source": [
    "i = start\n",
    "s = perf_counter()\n",
    "print(f\"Trajectory {i + 1}:\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 2.3 根据 Gaussian 分布产生随机的共轭动量场\n",
    "\n",
    "$$\n",
    "\\pi_\\mu(x)\\sim e^{-\\frac{1}{2}\\mathrm{Tr}\\left[\\pi_\\mu(x)\\pi_\\mu(x)\\right]}\\\\\n",
    "\\pi^\\dagger_\\mu(x)=\\pi_\\mu(x),\\mathrm{Tr}\\left[\\pi_\\mu(x)\\right]=0\n",
    "$$\n",
    "\n",
    "使用 Box–Muller 变换从均匀分布随机数产生高斯分布随机数。\n",
    "\n",
    "QUDA 实际产生 $i\\pi_\\mu(x)$ 作为动量场，可以通过 `saveMom()` 取出。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Creating Gaussian distributed Lie algebra field\n"
     ]
    }
   ],
   "source": [
    "hmc.gaussMom(i)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 2.4 计算 Hamiltonion（体系能量）\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathcal{H}(\\pi_\\mu,U_\\mu)=&\\frac{1}{2}\\sum_{x,\\mu}\\mathrm{Tr}\\left[\\pi_\\mu(x)\\pi_\\mu(x)\\right]\\\\\n",
    "&+\\frac{\\beta}{N_c}\\sum_{x,\\mu,\\nu\\gt\\mu}\\mathfrak{Re}\\mathrm{Tr}\\left[\\mathbb{I}_{N_c}-U_\\mu(x)U_\\nu(x+\\hat{\\mu})U^\\dagger_\\mu(x+\\hat{\\nu})U^\\dagger_\\nu(x)\\right]\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "`actionMom()` 用于计算动能项。\n",
    "\n",
    "`actionGauge()` 用于计算势能项。此时需要输入表示作用量具体组成的 `input_path` 与 `input_coeff`。\n",
    "\n",
    "此处动能和势能项可以添加任意常数，QUDA 实际使用下式计算总能量。\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathcal{H}(\\pi_\\mu,U_\\mu)=&\\frac{1}{2}\\sum_{x,\\mu}\\left\\{\\mathrm{Tr}\\left[\\pi_\\mu(x)\\pi_\\mu(x)\\right]-8\\right\\}\\\\\n",
    "&-\\frac{\\beta}{N_c}\\sum_{x,\\mu,\\nu\\gt\\mu}\\mathfrak{Re}\\mathrm{Tr}\\left[U_\\mu(x)U_\\nu(x+\\hat{\\mu})U^\\dagger_\\mu(x+\\hat{\\nu})U^\\dagger_\\nu(x)\\right]\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P_old = -4875878.399999999, K_old = 908.678913355529\n",
      "\n"
     ]
    }
   ],
   "source": [
    "kinetic_old = hmc.actionMom()\n",
    "potential_old = hmc.actionGauge(input_path, input_coeff)\n",
    "energy_old = kinetic_old + potential_old\n",
    "print(f\"P_old = {potential_old}, K_old = {kinetic_old}\\n\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "908.6789133553393\n",
      "-4875878.4\n"
     ]
    }
   ],
   "source": [
    "tmp_gauge = core.LatticeGauge(latt_info)\n",
    "tmp_mom = core.LatticeMom(latt_info)\n",
    "hmc.saveGauge(tmp_gauge)\n",
    "hmc.saveMom(tmp_mom)\n",
    "\n",
    "print(\n",
    "    (tmp_mom.data[:, :, :, :, :, :, :6] ** 2).sum()\n",
    "    + (tmp_mom.data[:, :, :, :, :, :, 6:] ** 2).sum() / 2\n",
    "    - 4 * 4 * latt_info.volume\n",
    ")\n",
    "print(-beta / Nc * hmc.plaquette() * Nc * 6 * latt_info.volume)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 2.5 分子动力学积分\n",
    "\n",
    "[Computer Physics Communications, Volume 151, Issue 3, 1 April 2003, Pages 272-314](https://www.sciencedirect.com/science/article/abs/pii/S0010465502007543)\n",
    "\n",
    "根据 Hamiltonian 正则方程，有分子动力学积分步骤\n",
    "$$\n",
    "\\pi'_\\mu(x)=\\pi_\\mu(x)+F_g(U_\\mu)\\delta t,\\;U'_\\mu(x)=e^{i\\pi_\\mu(x)\\delta t}U_\\mu(x)\n",
    "$$\n",
    "将其分为更新规范场 $U_\\mu(x)$ 的 $\\mathcal{A}$ 步骤和共轭动量场 $\\pi_\\mu(x)$ 的 $\\mathcal{B}$ 步骤。\n",
    "\n",
    "对于 Hamiltonian 体系的微正则系综，需要采用保辛结构积分保证总能量守恒。假设积分步长为 $\\Delta t$，一个 2 阶的保辛结构积分形式为\n",
    "$$\n",
    "\\mathcal{B}(\\Delta t/2)\\mathcal{A}(\\Delta t)\\mathcal{B}(\\Delta t/2)\n",
    "$$\n",
    "数值算法无法取到无穷小的 $\\Delta t$，导致积分完成后总能量出现一个小的偏离 $\\Delta\\mathcal{H}\\sim\\mathcal{O}(\\Delta t^3)$。通常我们会连续进行 $N$ 次 $\\Delta t$ 积分直到积分时间为 $\\tau$，那么总的偏离来到 $\\mathcal{O}(\\tau\\Delta t^2)=\\mathcal{O}(\\tau^3N^{-2})$。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "dt = t / n_steps\n",
    "for _ in range(n_steps):\n",
    "    hmc.updateMom(input_path2, input_coeff2, dt / 2)\n",
    "    hmc.updateGauge(dt)\n",
    "    hmc.updateMom(input_path2, input_coeff2, dt / 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 2.6 规范场重投影\n",
    "\n",
    "由于数值计算的舍入误差，规范场将会略微偏离 SU(3)，需要 `reunitGauge()` 重新投影。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [],
   "source": [
    "hmc.reunitGauge(1e-15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 2.7 计算新的 Hamiltonion（体系能量）\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\mathcal{H}(\\pi'_\\mu,U'_\\mu)=&\\frac{1}{2}\\sum_{x,\\mu}\\left\\{\\mathrm{Tr}\\left[\\pi'_\\mu(x)\\pi'_\\mu(x)\\right]-8\\right\\}\\\\\n",
    "&-\\frac{\\beta}{N_c}\\sum_{x,\\mu,\\nu\\gt\\mu}\\mathfrak{Re}\\mathrm{Tr}\\left[U'_\\mu(x)U'_\\nu(x+\\hat{\\mu})U^{\\prime\\dagger}_\\mu(x+\\hat{\\nu})U^{\\prime\\dagger}_\\nu(x)\\right]\n",
    "\\end{aligned}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "P = -4136377.3380672657, K = -729554.9059053268\n",
      "\n"
     ]
    }
   ],
   "source": [
    "kinetic = hmc.actionMom()\n",
    "potential = hmc.actionGauge(input_path, input_coeff)\n",
    "energy = kinetic + potential\n",
    "print(f\"P = {potential}, K = {kinetic}\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 2.8 判断是否接受新的规范场\n",
    "\n",
    "考虑 Markov 链的细致平衡条件\n",
    "$$\n",
    "P(U'_\\mu|U_\\mu)\\pi(U_\\mu)=P(U_\\mu|U'_\\mu)\\pi(U'_\\mu)\n",
    "$$\n",
    "\n",
    "将状态转移概率 $P(U'_\\mu|U_\\mu)$ 拆分为提议分布和接受概率 $g(U'_\\mu|U_\\mu)A(U'_\\mu|U_\\mu)$\n",
    "$$\n",
    "\\frac{A(U'_\\mu|U_\\mu)}{A(U_\\mu|U'_\\mu)}=\\frac{g(U_\\mu|U'_\\mu)\\pi(U'_\\mu)}{g(U'_\\mu|U_\\mu)\\pi(U_\\mu)}\n",
    "$$\n",
    "\n",
    "将 Hamiltonian 动能项记为 $\\mathcal{K}$，势能项记为 $\\mathcal{P}$，那么 $\\mathcal{H}=\\mathcal{K}+\\mathcal{P}$。\n",
    "\n",
    "对于分子动力学积分步骤，两个态互相转变分别需要共轭动量 $\\pi_\\mu$ 和 $-\\pi'_\\mu$，那么提议分布有比值\n",
    "$$\n",
    "\\frac{g(U_\\mu|U'_\\mu)}{g(U'_\\mu|U_\\mu)}=\\frac{e^{-\\mathcal{K}'}}{e^{-\\mathcal{K}}}\n",
    "$$\n",
    "平衡时，态密度 $\\pi(U_\\mu)=e^{-\\mathcal{P}}$\n",
    "$$\n",
    "\\frac{A(U'_\\mu|U_\\mu)}{A(U_\\mu|U'_\\mu)}=\\frac{e^{-\\mathcal{K}'}\\pi(U'_\\mu)}{e^{-\\mathcal{K}}\\pi(U_\\mu)}=e^{-\\mathcal{H}'+\\mathcal{H}}\n",
    "$$\n",
    "\n",
    "Metropolis 算法给出一种接受率选择\n",
    "$$\n",
    "A(U'_\\mu|U_\\mu)=\\min\\left[1,e^{-\\mathcal{H}'+\\mathcal{H}}\\right]\n",
    "$$\n",
    "\n",
    "`random()` 产生一个 $[0, 1)$ 区间内均匀分布的浮点数 $a$，通过 $a<A(\\pi'_\\mu,U'_\\mu;\\pi_\\mu,U_\\mu)$ 决定是否接受新的状态。当前步数低于 `warm` 时，永远接受新的状态。\n",
    "\n",
    "若接受，则通过 `saveGauge()` 保存当前规范场；若拒绝，则通过 `loadGauge()` 恢复旧的规范场。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Delta_P = 739501.0619327337, Delta_K = -730463.5848186824\n",
      "Delta_E = 9037.47711405158\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "\n"
     ]
    }
   ],
   "source": [
    "accept = random() < exp(energy_old - energy)\n",
    "if accept or i < warm:\n",
    "    hmc.saveGauge(gauge)\n",
    "else:\n",
    "    hmc.loadGauge(gauge)\n",
    "print(\n",
    "    f\"Delta_P = {potential - potential_old}, Delta_K = {kinetic - kinetic_old}\\n\"\n",
    "    f\"Delta_E = {energy - energy_old}\\n\"\n",
    "    f\"acceptance rate = {min(1, exp(energy_old - energy))*100:.2f}%\\n\"\n",
    "    f\"accept? {accept or i < warm}\\n\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 2.9 第 1 步结束，输出统计信息\n",
    "\n",
    "根据 `save` 参数决定是否保存当前规范场为组态文件。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Trajectory 1:\n",
      "plaquette = 0.8483348022106675\n",
      "P_old = -4875878.399999999, K_old = 908.678913355529\n",
      "P = -4136377.3380672657, K = -729554.9059053268\n",
      "Delta_P = 739501.0619327337, Delta_K = -730463.5848186824\n",
      "Delta_E = 9037.47711405158\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 30.937 secs\n",
      "\n"
     ]
    }
   ],
   "source": [
    "plaquette.append(hmc.plaquette())\n",
    "print(\n",
    "    f\"Trajectory {i + 1}:\\n\"\n",
    "    f\"plaquette = {plaquette[-1]}\\n\"\n",
    "    f\"P_old = {potential_old}, K_old = {kinetic_old}\\n\"\n",
    "    f\"P = {potential}, K = {kinetic}\\n\"\n",
    "    f\"Delta_P = {potential - potential_old}, Delta_K = {kinetic - kinetic_old}\\n\"\n",
    "    f\"Delta_E = {energy - energy_old}\\n\"\n",
    "    f\"acceptance rate = {min(1, exp(energy_old - energy))*100:.2f}%\\n\"\n",
    "    f\"accept? {accept or i < warm}\\n\"\n",
    "    f\"HMC time = {perf_counter() - s:.3f} secs\\n\"\n",
    ")\n",
    "\n",
    "if (i + 1) % save == 0:\n",
    "    io.writeNPYGauge(f\"./DATA/cfg/cfg_{i + 1}.npy\", gauge)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 2.10 循环 2.2～2.9\n",
    "\n",
    "循环从 `start` 开始，至 `stop` 结束。其中第 1 步在循环外完成，计数器从 `start + 1` 开始。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 2:\n",
      "plaquette = 0.7825983901313961\n",
      "P_old = -4136377.3380672657, K_old = -1005.1336467889819\n",
      "P = -3815854.5863164477, K = -318504.96315335215\n",
      "Delta_P = 320522.751750818, Delta_K = -317499.82950656314\n",
      "Delta_E = 3022.9222442549653\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.138 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 3:\n",
      "plaquette = 0.7363945870715864\n",
      "P_old = -3815854.5863164477, K_old = -344.00956808673914\n",
      "P = -3590570.4609792675, K = -223697.42204911186\n",
      "Delta_P = 225284.12533718022, Delta_K = -223353.41248102512\n",
      "Delta_E = 1930.7128561553545\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 4:\n",
      "plaquette = 0.7019395394857717\n",
      "P_old = -3590570.4609792675, K_old = 4039.1828641460575\n",
      "P = -3422571.838684622, K = -162616.4618685979\n",
      "Delta_P = 167998.6222946453, Delta_K = -166655.64473274394\n",
      "Delta_E = 1342.9775619013235\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 5:\n",
      "plaquette = 0.6767887444930591\n",
      "P_old = -3422571.838684622, K_old = -725.7948669571351\n",
      "P = -3299939.6206368264, K = -122424.07304089903\n",
      "Delta_P = 122632.21804779582, Delta_K = -121698.2781739419\n",
      "Delta_E = 933.939873854164\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 6:\n",
      "plaquette = 0.6583577662110424\n",
      "P_old = -3299939.6206368264, K_old = -825.3822401503437\n",
      "P = -3210072.411740671, K = -90023.26197021533\n",
      "Delta_P = 89867.20889615547, Delta_K = -89197.87973006499\n",
      "Delta_E = 669.3291660901159\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 7:\n",
      "plaquette = 0.6446671479394226\n",
      "P_old = -3210072.411740671, K_old = -633.59363864513\n",
      "P = -3143318.6218274357, K = -66903.0827277755\n",
      "Delta_P = 66753.78991323523, Delta_K = -66269.48908913037\n",
      "Delta_E = 484.30082410480827\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 8:\n",
      "plaquette = 0.6350341432927649\n",
      "P_old = -3143318.6218274357, K_old = -242.57249829810826\n",
      "P = -3096349.2625436983, K = -46874.860087966\n",
      "Delta_P = 46969.35928373737, Delta_K = -46632.28758966789\n",
      "Delta_E = 337.071694069542\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 9:\n",
      "plaquette = 0.62804789123982\n",
      "P_old = -3096349.2625436983, K_old = 1010.4361033365166\n",
      "P = -3062285.1470617875, K = -32819.24685225683\n",
      "Delta_P = 34064.11548191076, Delta_K = -33829.682955593344\n",
      "Delta_E = 234.43252631742507\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.130 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 10:\n",
      "plaquette = 0.6229850354817332\n",
      "P_old = -3062285.1470617875, K_old = -562.3321330886077\n",
      "P = -3037599.2780286167, K = -25063.34505319116\n",
      "Delta_P = 24685.869033170864, Delta_K = -24501.012920102552\n",
      "Delta_E = 184.85611306829378\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 11:\n",
      "plaquette = 0.6191817307243552\n",
      "P_old = -3037599.2780286167, K_old = -1453.8689941094674\n",
      "P = -3019054.8265135, K = -19868.943893135383\n",
      "Delta_P = 18544.45151511673, Delta_K = -18415.074899025916\n",
      "Delta_E = 129.37661609053612\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 12:\n",
      "plaquette = 0.6165636670992021\n",
      "P_old = -3019054.8265135, K_old = -467.6170520559517\n",
      "P = -3006289.46663379, K = -13132.16685043215\n",
      "Delta_P = 12765.35987970978, Delta_K = -12664.549798376198\n",
      "Delta_E = 100.81008133338764\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.130 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 13:\n",
      "plaquette = 0.6145003041001613\n",
      "P_old = -3006289.46663379, K_old = 606.3785938293045\n",
      "P = -2996228.7595554083, K = -9386.93706136494\n",
      "Delta_P = 10060.707078381907, Delta_K = -9993.315655194245\n",
      "Delta_E = 67.39142318768427\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 14:\n",
      "plaquette = 0.6132765064617612\n",
      "P_old = -2996228.7595554083, K_old = 1003.9146727464869\n",
      "P = -2990261.6710843625, K = -4896.427513729735\n",
      "Delta_P = 5967.088471045718, Delta_K = -5900.342186476222\n",
      "Delta_E = 66.74628456961364\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 15:\n",
      "plaquette = 0.6128130873213314\n",
      "P_old = -2990261.6710843625, K_old = -1025.2754457459332\n",
      "P = -2988002.095707394, K = -3253.684225620941\n",
      "Delta_P = 2259.575376968365, Delta_K = -2228.408779875008\n",
      "Delta_E = 31.166597092989832\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 16:\n",
      "plaquette = 0.6124003181723601\n",
      "P_old = -2988002.095707394, K_old = 57.8440681975419\n",
      "P = -2985989.4835297386, K = -1936.0605669063575\n",
      "Delta_P = 2012.6121776555665, Delta_K = -1993.9046351038994\n",
      "Delta_E = 18.7075425516814\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.130 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 17:\n",
      "plaquette = 0.6120091618202393\n",
      "P_old = -2985989.4835297386, K_old = 136.27114785181948\n",
      "P = -2984082.2527214102, K = -1737.4487564206515\n",
      "Delta_P = 1907.2308083283715, Delta_K = -1873.719904272471\n",
      "Delta_E = 33.51090405602008\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 18:\n",
      "plaquette = 0.6116508912202814\n",
      "P_old = -2984082.2527214102, K_old = 875.1615478352022\n",
      "P = -2982335.36884172, K = -845.4448556402926\n",
      "Delta_P = 1746.883879690431, Delta_K = -1720.606403475495\n",
      "Delta_E = 26.277476215269417\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.130 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 19:\n",
      "plaquette = 0.6113557669692254\n",
      "P_old = -2982335.36884172, K_old = -1060.8341325009164\n",
      "P = -2980896.3788806796, K = -2471.874145424583\n",
      "Delta_P = 1438.9899610402063, Delta_K = -1411.0400129236664\n",
      "Delta_E = 27.949948116671294\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 20:\n",
      "plaquette = 0.6110346513647702\n",
      "P_old = -2980896.3788806796, K_old = -2743.2553047591036\n",
      "P = -2979330.658241014, K = -4288.431690802512\n",
      "Delta_P = 1565.720639665611, Delta_K = -1545.1763860434085\n",
      "Delta_E = 20.544253622181714\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.130 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 21:\n",
      "plaquette = 0.6108110894443504\n",
      "P_old = -2979330.658241014, K_old = 466.9325671204556\n",
      "P = -2978240.597502176, K = -612.8740309586879\n",
      "Delta_P = 1090.060738837812, Delta_K = -1079.8065980791434\n",
      "Delta_E = 10.254140758421272\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 22:\n",
      "plaquette = 0.610660370995803\n",
      "P_old = -2978240.597502176, K_old = -475.82221571335253\n",
      "P = -2977505.712674422, K = -1189.984307959947\n",
      "Delta_P = 734.8848277539946, Delta_K = -714.1620922465945\n",
      "Delta_E = 20.722735507413745\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.130 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 23:\n",
      "plaquette = 0.6109364169053637\n",
      "P_old = -2977505.712674422, K_old = -983.0152892393011\n",
      "P = -2978851.6789622577, K = 367.2959765698607\n",
      "Delta_P = -1345.9662878355011, Delta_K = 1350.3112658091618\n",
      "Delta_E = 4.344977973494679\n",
      "acceptance rate = 1.30%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 24:\n",
      "plaquette = 0.6105514112550297\n",
      "P_old = -2978851.6789622577, K_old = -632.6174889643116\n",
      "P = -2976974.438227915, K = -2479.9727889822416\n",
      "Delta_P = 1877.2407343424857, Delta_K = -1847.35530001793\n",
      "Delta_E = 29.885434324853122\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.130 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 25:\n",
      "plaquette = 0.609813685960289\n",
      "P_old = -2976974.438227915, K_old = -24.867148569352537\n",
      "P = -2973377.3793981574, K = -3580.7366034261086\n",
      "Delta_P = 3597.0588297578506, Delta_K = -3555.869454856756\n",
      "Delta_E = 41.18937490088865\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 26:\n",
      "plaquette = 0.609585541445256\n",
      "P_old = -2973377.3793981574, K_old = 1328.263784522278\n",
      "P = -2972264.9744852283, K = 232.2308480718409\n",
      "Delta_P = 1112.4049129290506, Delta_K = -1096.032936450437\n",
      "Delta_E = 16.37197647849098\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.130 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 27:\n",
      "plaquette = 0.6098066697698308\n",
      "P_old = -2972264.9744852283, K_old = -1848.521658016445\n",
      "P = -2973343.1693066508, K = -770.1852092597327\n",
      "Delta_P = -1078.194821422454, Delta_K = 1078.3364487567123\n",
      "Delta_E = 0.14162733405828476\n",
      "acceptance rate = 86.79%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 28:\n",
      "plaquette = 0.6101012398795089\n",
      "P_old = -2973343.1693066508, K_old = -171.80966122414787\n",
      "P = -2974779.4573417166, K = 1269.8691323512799\n",
      "Delta_P = -1436.288035065867, Delta_K = 1441.6787935754278\n",
      "Delta_E = 5.390758509282023\n",
      "acceptance rate = 0.46%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 29:\n",
      "plaquette = 0.6103196899212264\n",
      "P_old = -2974779.4573417166, K_old = -504.81587146901495\n",
      "P = -2975844.5931816064, K = 570.9698866649459\n",
      "Delta_P = -1065.1358398897573, Delta_K = 1075.7857581339608\n",
      "Delta_E = 10.649918244220316\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 30:\n",
      "plaquette = 0.6101288468627317\n",
      "P_old = -2975844.5931816064, K_old = -194.69796496464596\n",
      "P = -2974914.0656349007, K = -1103.7897104868357\n",
      "Delta_P = 930.5275467056781, Delta_K = -909.0917455221897\n",
      "Delta_E = 21.43580118333921\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.130 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 31:\n",
      "plaquette = 0.6101545556054748\n",
      "P_old = -2974914.0656349007, K_old = 1481.3596436047883\n",
      "P = -2975039.418338334, K = 1617.2964212606175\n",
      "Delta_P = -125.35270343348384, Delta_K = 135.93677765582925\n",
      "Delta_E = 10.58407422248274\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 32:\n",
      "plaquette = 0.6100866364913988\n",
      "P_old = -2975039.418338334, K_old = 1216.487128687495\n",
      "P = -2974708.252997063, K = 895.0703080908856\n",
      "Delta_P = 331.1653412710875, Delta_K = -321.41682059660934\n",
      "Delta_E = 9.748520674183965\n",
      "acceptance rate = 0.01%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 33:\n",
      "plaquette = 0.6105557448053432\n",
      "P_old = -2974708.252997063, K_old = -1019.9748483753394\n",
      "P = -2976995.5680922857, K = 1261.222865911259\n",
      "Delta_P = -2287.315095222555, Delta_K = 2281.1977142865985\n",
      "Delta_E = -6.117380936164409\n",
      "acceptance rate = 100.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 34:\n",
      "plaquette = 0.6102470367988171\n",
      "P_old = -2976995.5680922857, K_old = 1263.963581728221\n",
      "P = -2975490.345391358, K = -220.3427218478951\n",
      "Delta_P = 1505.2227009278722, Delta_K = -1484.306303576116\n",
      "Delta_E = 20.916397351771593\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.130 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 35:\n",
      "plaquette = 0.6104908924886018\n",
      "P_old = -2975490.345391358, K_old = -3911.9134634345924\n",
      "P = -2976679.356081896, K = -2716.323844917762\n",
      "Delta_P = -1189.0106905382127, Delta_K = 1195.5896185168303\n",
      "Delta_E = 6.578927978873253\n",
      "acceptance rate = 0.14%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 36:\n",
      "plaquette = 0.6107058915628683\n",
      "P_old = -2976679.356081896, K_old = -1870.4696299301013\n",
      "P = -2977727.6654241313, K = -821.9457590494694\n",
      "Delta_P = -1048.3093422353268, Delta_K = 1048.5238708806319\n",
      "Delta_E = 0.21452864538878202\n",
      "acceptance rate = 80.69%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 37:\n",
      "plaquette = 0.6106622028558144\n",
      "P_old = -2977727.6654241313, K_old = -2389.512655396878\n",
      "P = -2977514.6446010843, K = -2590.6771426120677\n",
      "Delta_P = 213.02082304703072, Delta_K = -201.16448721518964\n",
      "Delta_E = 11.856335831806064\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 38:\n",
      "plaquette = 0.6109228282071942\n",
      "P_old = -2977514.6446010843, K_old = -2094.6554130516547\n",
      "P = -2978785.422122369, K = -802.8457227503608\n",
      "Delta_P = -1270.7775212847628, Delta_K = 1291.8096903012938\n",
      "Delta_E = 21.032169016543776\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.130 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 39:\n",
      "plaquette = 0.610732159098142\n",
      "P_old = -2978785.422122369, K_old = 1240.922244940269\n",
      "P = -2977855.742731994, K = 319.98669569700894\n",
      "Delta_P = 929.6793903750367, Delta_K = -920.93554924326\n",
      "Delta_E = 8.743841131683439\n",
      "acceptance rate = 0.02%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 40:\n",
      "plaquette = 0.6105738016667808\n",
      "P_old = -2977855.742731994, K_old = 150.06769597463045\n",
      "P = -2977083.6111529414, K = -598.0428303248473\n",
      "Delta_P = 772.1315790526569, Delta_K = -748.1105262994777\n",
      "Delta_E = 24.021052753087133\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.130 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 41:\n",
      "plaquette = 0.6103211439246816\n",
      "P_old = -2977083.6111529414, K_old = 2376.5970747868764\n",
      "P = -2975851.6827256465, K = 1155.8971510467204\n",
      "Delta_P = 1231.928427294828, Delta_K = -1220.699923740156\n",
      "Delta_E = 11.228503555059433\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 42:\n",
      "plaquette = 0.6104024094157772\n",
      "P_old = -2975851.6827256465, K_old = -1108.4327265563757\n",
      "P = -2976247.923378345, K = -698.1330291073516\n",
      "Delta_P = -396.2406526985578, Delta_K = 410.29969744902405\n",
      "Delta_E = 14.059044750407338\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 43:\n",
      "plaquette = 0.609986299523507\n",
      "P_old = -2976247.923378345, K_old = 1576.7131331970113\n",
      "P = -2974219.0221425984, K = -426.25023802774666\n",
      "Delta_P = 2028.9012357466854, Delta_K = -2002.963371224758\n",
      "Delta_E = 25.93786452198401\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 44:\n",
      "plaquette = 0.61007490620552\n",
      "P_old = -2974219.0221425984, K_old = 833.8018127669429\n",
      "P = -2974651.057549521, K = 1272.5534207413725\n",
      "Delta_P = -432.0354069224559, Delta_K = 438.7516079744296\n",
      "Delta_E = 6.716201052069664\n",
      "acceptance rate = 0.12%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 45:\n",
      "plaquette = 0.6104721277796405\n",
      "P_old = -2974651.057549521, K_old = -57.570348864050175\n",
      "P = -2976587.8616427896, K = 1886.911268886673\n",
      "Delta_P = -1936.8040932687, Delta_K = 1944.4816177507232\n",
      "Delta_E = 7.677524481900036\n",
      "acceptance rate = 0.05%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 46:\n",
      "plaquette = 0.6103181166547427\n",
      "P_old = -2976587.8616427896, K_old = -93.81428742627418\n",
      "P = -2975836.9221255397, K = -821.595441087667\n",
      "Delta_P = 750.9395172498189, Delta_K = -727.7811536613929\n",
      "Delta_E = 23.158363588619977\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.130 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 47:\n",
      "plaquette = 0.6107740659455987\n",
      "P_old = -2975836.9221255397, K_old = -4678.35580895125\n",
      "P = -2978060.0754243196, K = -2450.724697636044\n",
      "Delta_P = -2223.1532987798564, Delta_K = 2227.631111315206\n",
      "Delta_E = 4.477812535595149\n",
      "acceptance rate = 1.14%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 48:\n",
      "plaquette = 0.6108782462865994\n",
      "P_old = -2978060.0754243196, K_old = -1409.8212072226806\n",
      "P = -2978568.04609871, K = -891.8321965954347\n",
      "Delta_P = -507.97067439043894, Delta_K = 517.989010627246\n",
      "Delta_E = 10.018336236476898\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 49:\n",
      "plaquette = 0.6104048828472793\n",
      "P_old = -2978568.04609871, K_old = 515.5573136203088\n",
      "P = -2976259.98352958, K = -1752.9866623447929\n",
      "Delta_P = 2308.062569130212, Delta_K = -2268.5439759651017\n",
      "Delta_E = 39.51859316509217\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n",
      "Creating Gaussian distributed Lie algebra field\n",
      "Trajectory 50:\n",
      "plaquette = 0.6100999421041797\n",
      "P_old = -2976259.98352958, K_old = 47.930575923478784\n",
      "P = -2974773.129547021, K = -1414.8321090645973\n",
      "Delta_P = 1486.8539825589396, Delta_K = -1462.762684988076\n",
      "Delta_E = 24.091297570616007\n",
      "acceptance rate = 0.00%\n",
      "accept? True\n",
      "HMC time = 0.131 secs\n",
      "\n"
     ]
    },
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: './DATA/cfg/cfg_50.npy'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[16], line 42\u001b[0m\n\u001b[1;32m     29\u001b[0m \u001b[38;5;28mprint\u001b[39m(\n\u001b[1;32m     30\u001b[0m     \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mTrajectory \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mi\u001b[38;5;250m \u001b[39m\u001b[38;5;241m+\u001b[39m\u001b[38;5;250m \u001b[39m\u001b[38;5;241m1\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m:\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m     31\u001b[0m     \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mplaquette = \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mplaquette[\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m]\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m   (...)\u001b[0m\n\u001b[1;32m     38\u001b[0m     \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mHMC time = \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mperf_counter()\u001b[38;5;250m \u001b[39m\u001b[38;5;241m-\u001b[39m\u001b[38;5;250m \u001b[39ms\u001b[38;5;132;01m:\u001b[39;00m\u001b[38;5;124m.3f\u001b[39m\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m secs\u001b[39m\u001b[38;5;130;01m\\n\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m     39\u001b[0m )\n\u001b[1;32m     41\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m (i \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39m) \u001b[38;5;241m%\u001b[39m save \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[0;32m---> 42\u001b[0m     \u001b[43mio\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwriteNPYGauge\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;124;43mf\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43m./DATA/cfg/cfg_\u001b[39;49m\u001b[38;5;132;43;01m{\u001b[39;49;00m\u001b[43mi\u001b[49m\u001b[38;5;250;43m \u001b[39;49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[38;5;250;43m \u001b[39;49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[38;5;132;43;01m}\u001b[39;49;00m\u001b[38;5;124;43m.npy\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mgauge\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/.venv/lib/python3.11/site-packages/pyquda_utils/io/__init__.py:205\u001b[0m, in \u001b[0;36mwriteNPYGauge\u001b[0;34m(filename, gauge)\u001b[0m\n\u001b[1;32m    202\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01m.\u001b[39;00m\u001b[38;5;21;01mnpy\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m writeGauge \u001b[38;5;28;01mas\u001b[39;00m write\n\u001b[1;32m    204\u001b[0m filename \u001b[38;5;241m=\u001b[39m filename \u001b[38;5;28;01mif\u001b[39;00m filename\u001b[38;5;241m.\u001b[39mendswith(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m.npy\u001b[39m\u001b[38;5;124m\"\u001b[39m) \u001b[38;5;28;01melse\u001b[39;00m filename \u001b[38;5;241m+\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m.npy\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m--> 205\u001b[0m \u001b[43mwrite\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfilename\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mgauge\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlatt_info\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mglobal_size\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mgauge\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mlexico\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m~/.venv/lib/python3.11/site-packages/pyquda_utils/io/npy.py:36\u001b[0m, in \u001b[0;36mwriteGauge\u001b[0;34m(filename, latt_size, gauge)\u001b[0m\n\u001b[1;32m     34\u001b[0m Lx, Ly, Lz, Lt \u001b[38;5;241m=\u001b[39m getSublatticeSize(latt_size)\n\u001b[1;32m     35\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m getMPIRank() \u001b[38;5;241m==\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[0;32m---> 36\u001b[0m     \u001b[38;5;28;01mwith\u001b[39;00m \u001b[38;5;28;43mopen\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mfilename\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mwb\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mas\u001b[39;00m f:\n\u001b[1;32m     37\u001b[0m         write_array_header_1_0(\n\u001b[1;32m     38\u001b[0m             f, {\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mshape\u001b[39m\u001b[38;5;124m\"\u001b[39m: (Nd, \u001b[38;5;241m*\u001b[39mlatt_size[::\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m], Nc, Nc), \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mfortran_order\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;28;01mFalse\u001b[39;00m, \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mdescr\u001b[39m\u001b[38;5;124m\"\u001b[39m: \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m<c16\u001b[39m\u001b[38;5;124m\"\u001b[39m}\n\u001b[1;32m     39\u001b[0m         )\n\u001b[1;32m     40\u001b[0m getMPIComm()\u001b[38;5;241m.\u001b[39mBarrier()\n",
      "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: './DATA/cfg/cfg_50.npy'"
     ]
    }
   ],
   "source": [
    "for i in range(start + 1, stop):\n",
    "    s = perf_counter()\n",
    "\n",
    "    hmc.gaussMom(i)\n",
    "\n",
    "    kinetic_old = hmc.actionMom()\n",
    "    potential_old = hmc.actionGauge(input_path, input_coeff)\n",
    "    energy_old = kinetic_old + potential_old\n",
    "\n",
    "    dt = t / n_steps\n",
    "    for _ in range(n_steps):\n",
    "        hmc.updateMom(input_path2, input_coeff2, dt / 2)\n",
    "        hmc.updateGauge(dt)\n",
    "        hmc.updateMom(input_path2, input_coeff2, dt / 2)\n",
    "\n",
    "    hmc.reunitGauge(1e-15)\n",
    "\n",
    "    kinetic = hmc.actionMom()\n",
    "    potential = hmc.actionGauge(input_path, input_coeff)\n",
    "    energy = kinetic + potential\n",
    "\n",
    "    accept = random() < exp(energy_old - energy)\n",
    "    if accept or i < warm:\n",
    "        hmc.saveGauge(gauge)\n",
    "    else:\n",
    "        hmc.loadGauge(gauge)\n",
    "\n",
    "    plaquette.append(hmc.plaquette())\n",
    "    print(\n",
    "        f\"Trajectory {i + 1}:\\n\"\n",
    "        f\"plaquette = {plaquette[-1]}\\n\"\n",
    "        f\"P_old = {potential_old}, K_old = {kinetic_old}\\n\"\n",
    "        f\"P = {potential}, K = {kinetic}\\n\"\n",
    "        f\"Delta_P = {potential - potential_old}, Delta_K = {kinetic - kinetic_old}\\n\"\n",
    "        f\"Delta_E = {energy - energy_old}\\n\"\n",
    "        f\"acceptance rate = {min(1, exp(energy_old - energy))*100:.2f}%\\n\"\n",
    "        f\"accept? {accept or i < warm}\\n\"\n",
    "        f\"HMC time = {perf_counter() - s:.3f} secs\\n\"\n",
    "    )\n",
    "\n",
    "    if (i + 1) % save == 0:\n",
    "        io.writeNPYGauge(f\"./DATA/cfg/cfg_{i + 1}.npy\", gauge)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### 2.11 画图"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.6, 0.62)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAGiCAYAAADnfswJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABOiUlEQVR4nO3de1xUZf4H8M/MwAyiXFScYUAQL6GYgoYyTm5ZibfsYrm7tlmav1YLRyWoVtldJbuALbuulWymm+nuWrq6maaoGZpWkhcU76CIiBcGRIRBFAZmnt8fyNQEKkPACOfzfr3O68Wc85znfM+Rmg/nPOccmRBCgIiIiEgC5M4ugIiIiKilMPgQERGRZDD4EBERkWQw+BAREZFkMPgQERGRZDD4EBERkWQw+BAREZFkMPgQERGRZDD4EBERkWQw+BAREZFkNCr4JCcnIygoCG5ubtDpdNi3b99t25eUlMBgMECr1UKlUiE4OBgpKSm25YmJiRg8eDA8PDygVqsxbtw4ZGVl2fVRUVEBg8GAzp07o0OHDhg/fjwKCgrs2uTl5WHs2LFwd3eHWq3G66+/jurq6sbsIhEREbVBDgefNWvWIDY2FvHx8Th48CDCwsIwatQoFBYW1tvebDZjxIgRyM3Nxbp165CVlYVly5bB39/f1mbXrl0wGAz44YcfsH37dlRVVWHkyJEoLy+3tYmJicGXX36JtWvXYteuXbh06RKefvpp23KLxYKxY8fCbDZjz549WLlyJVasWIF58+Y5uotERETUVgkHRURECIPBYPtssViEn5+fSExMrLf9hx9+KHr06CHMZnODt1FYWCgAiF27dgkhhCgpKRGurq5i7dq1tjYnT54UAERaWpoQQoiUlBQhl8uF0Wi027anp6eorKx0aB+JiIiobXJxJCSZzWakp6cjLi7ONk8ulyMyMhJpaWn1rrNx40bo9XoYDAZs2LABXbp0wbPPPovZs2dDoVDUu05paSkAoFOnTgCA9PR0VFVVITIy0tamT58+CAwMRFpaGoYMGYK0tDT0798fGo3G1mbUqFGIiorC8ePHMXDgwDrbqaysRGVlpe2z1WpFcXExOnfuDJlM5sCRISIiImcRQqCsrAx+fn6Qy29/Mcuh4FNUVASLxWIXLgBAo9EgMzOz3nVycnKwY8cOTJw4ESkpKcjOzsb06dNRVVWF+Pj4Ou2tViteeeUVDB06FP369QMAGI1GKJVKeHt719mu0Wi0tamvrtpl9UlMTMT8+fPvvONERER01zt//jy6du162zYOBZ/GsFqtUKvVWLp0KRQKBcLDw3Hx4kUkJSXVG3wMBgOOHTuG7777rrlLQ1xcHGJjY22fS0tLERgYiPPnz8PT07PZt08t64F3d+Dq9Sqse1mPPlr++xIRtRUmkwkBAQHw8PC4Y1uHgo+Pjw8UCkWdu6kKCgrg6+tb7zparRaurq52l7VCQkJgNBphNpuhVCpt82fMmIFNmzZh9+7ddonN19cXZrMZJSUldmd9frpdX1/fOneX1dZ5q9pUKhVUKlWd+Z6engw+bVCQ1gelF0pxtdqF/75ERG1QQ4apOHRXl1KpRHh4OFJTU23zrFYrUlNTodfr611n6NChyM7OhtVqtc07deoUtFqtLfQIITBjxgysX78eO3bsQPfu3e36CA8Ph6urq912s7KykJeXZ9uuXq/H0aNH7e4u2759Ozw9PdG3b19HdpPaqK6d3AEA54uvO7kSIiJyFodvZ4+NjcWyZcuwcuVKnDx5ElFRUSgvL8eUKVMAAJMmTbIb/BwVFYXi4mJER0fj1KlT2Lx5MxISEmAwGGxtDAYD/vOf/+DTTz+Fh4cHjEYjjEYjbty4AQDw8vLCiy++iNjYWOzcuRPp6emYMmUK9Ho9hgwZAgAYOXIk+vbti+effx6HDx/Gtm3b8Oc//xkGg6HeszokPQEda4LPhas3nFwJERE5i8NjfCZMmIDLly9j3rx5MBqNGDBgALZu3WobSJyXl2c3ojogIADbtm1DTEwMQkND4e/vj+joaMyePdvW5sMPPwQAPPTQQ3bb+uSTT/DCCy8AAP7+979DLpdj/PjxqKysxKhRo/CPf/zD1lahUGDTpk2IioqCXq9H+/btMXnyZLz55puO7iK1UQGd2gHgGR8iIimTCSGEs4u4W5hMJnh5eaG0tJRjQNqgXacuY/LyfQjWdMBXMcOcXQ4RETURR76/+a4ukoyAjjVnfC5cvQHmfSIiaWLwIcnw79gOMhlw3WxBcbnZ2eUQEZETMPiQZKhcFNB4uAEAznOAMxGRJDH4kKRwgDMRkbQx+JCk1N7Sfv4qgw8RkRQx+JCkdO1Ye8aHl7qIiKSIwYckpfbpzRd4xoeISJIYfEhSai91fXu6yMmVEBGRMzD4kKT0UncAgJu3tVc7uRoiImppDD4kKV08VOjioYIQQKaxzNnlEBFRC2PwIckJ0dY8zvxkvsnJlRARUUtj8CHJ6Xsz+Jy4xOBDRCQ1DD4kOX39bgYfnvEhIpIcBh+SnNozPlnGMlisfFkpEZGUMPiQ5HT3aQ83Vzmumy04d6Xc2eUQEVELYvAhyVHIZejty8tdRERSxOBDksQBzkRE0sTgQ5LUV+sBgLe0ExFJDYMPSRLv7CIikiYGH5Kk3r6ekMmAAlMlrlyrdHY5RETUQhh8SJI6qFwQ1Lk9AOBkPl9dQUQkFQw+JFkhN8f5nMgvdXIlRETUUhh8SLJ4ZxcRkfQw+JBkcYAzEZH0MPiQZNW+pf3M5XJUVFmcXA0REbUEBh+SLF9PN3R0d4XFKnC64JqzyyEiohbA4EOSJZPJfnK5iwOciYikgMGHJK12gDNvaScikgYGH5K0EN7ZRUQkKQw+JGk/vbPLahVOroaIiJobgw9JWs8uHaBUyHGtshoXrt5wdjlERNTMGHxI0lwVctyj6QCAz/MhIpICBh+SPNsTnBl8iIjaPAYfkjzbOB8OcCYiavMYfEjyfrylncGHiKita1TwSU5ORlBQENzc3KDT6bBv377bti8pKYHBYIBWq4VKpUJwcDBSUlJsy3fv3o3HH38cfn5+kMlk+OKLL+r0IZPJ6p2SkpJsbYKCguosX7BgQWN2kSSkz83gc7HkBkqvVzm5GiIiak4OB581a9YgNjYW8fHxOHjwIMLCwjBq1CgUFhbW295sNmPEiBHIzc3FunXrkJWVhWXLlsHf39/Wpry8HGFhYUhOTr7ldvPz8+2m5cuXQyaTYfz48Xbt3nzzTbt2M2fOdHQXSWK82rmia8d2ADjOh4iorXNxdIWFCxdi6tSpmDJlCgBgyZIl2Lx5M5YvX445c+bUab98+XIUFxdjz549cHV1BVBzZuanxowZgzFjxtx2u76+vnafN2zYgIcffhg9evSwm+/h4VGnLdGd9NV64sLVGziRb4K+Z2dnl0NERM3EoTM+ZrMZ6enpiIyM/LEDuRyRkZFIS0urd52NGzdCr9fDYDBAo9GgX79+SEhIgMXS+LdhFxQUYPPmzXjxxRfrLFuwYAE6d+6MgQMHIikpCdXV1bfsp7KyEiaTyW4iaQrhOB8iIklw6IxPUVERLBYLNBqN3XyNRoPMzMx618nJycGOHTswceJEpKSkIDs7G9OnT0dVVRXi4+MbVfTKlSvh4eGBp59+2m7+rFmzcN9996FTp07Ys2cP4uLikJ+fj4ULF9bbT2JiIubPn9+oGqht4Z1dRETS4PClLkdZrVao1WosXboUCoUC4eHhuHjxIpKSkhodfJYvX46JEyfCzc3Nbn5sbKzt59DQUCiVSrz00ktITEyESqWq009cXJzdOiaTCQEBAY2qiVq32ju7TheWwVxthdKFNzwSEbVFDgUfHx8fKBQKFBQU2M0vKCi45bgarVYLV1dXKBQK27yQkBAYjUaYzWYolUqHCv7222+RlZWFNWvW3LGtTqdDdXU1cnNz0bt37zrLVSpVvYGIpKdrx3bwcHNBWUU1zly+Zrv0RUREbYtDf9YqlUqEh4cjNTXVNs9qtSI1NRV6vb7edYYOHYrs7GxYrVbbvFOnTkGr1TocegDg448/Rnh4OMLCwu7YNiMjA3K5HGq12uHtkLTIZDK+qZ2ISAIcPp8fGxuLZcuWYeXKlTh58iSioqJQXl5uu8tr0qRJiIuLs7WPiopCcXExoqOjcerUKWzevBkJCQkwGAy2NteuXUNGRgYyMjIAAGfPnkVGRgby8vLstm0ymbB27Vr8/ve/r1NXWloaFi1ahMOHDyMnJwerVq1CTEwMnnvuOXTs2NHR3SQJ4qsriIjaPofH+EyYMAGXL1/GvHnzYDQaMWDAAGzdutU24DkvLw9y+Y95KiAgANu2bUNMTAxCQ0Ph7++P6OhozJ4929bmwIEDePjhh22fa8fdTJ48GStWrLDNX716NYQQ+N3vflenLpVKhdWrV+ONN95AZWUlunfvjpiYGLsxPES3wwHORERtn0wIIZxdxN3CZDLBy8sLpaWl8PTkGA+pOXaxFI998B283V1xaO4IyGQyZ5dEREQN4Mj3N29dIbqpl7oDXOQylFyvQn5phbPLISKiZsDgQ3STm6sCvdQdAPByFxFRW8XgQ/QTfFM7EVHbxuBD9BMhvLOLiKhNY/Ah+gnbnV0MPkREbRKDD9FP1J7xOXflOsoqqpxcDRERNTUGH6Kf6NReCV/PmnfAZRnLnFwNERE1NQYfop/h5S4ioraLwYfoZ/rynV1ERG0Wgw/Rz9SO81m9/7yTKyEioqbG4EP0M/38a4KPq0KGiiqLk6shIqKmxOBD9DOBndzRxUOFKovA4fMlzi6HiIiaEIMP0c/IZDJEBHUCABw4d9XJ1RARUVNi8CGqx6CgjgCAfWeLnVwJERE1JQYfonoMvnnG5+C5q7BYhZOrISKipsLgQ1SPEK0nOqhcUFZZjUwjb2snImorGHyI6qGQy3Bft5rLXft5uYuIqM1g8CG6hcG1wYcDnImI2gwGH6JbGNy9ZpzP/rPFEILjfIiI2gIGH6JbGBDgDVeFDIVllThffMPZ5RARURNg8CG6BTdXBfr7ewEA9uVynA8RUVvA4EN0G7WXuw4w+BARtQkMPkS3MbhbTfDhGR8ioraBwYfoNmqf4JxzuRxXrlU6uRoiIvqlGHyIbsPbXYlgTQcAwP5c3tZORNTaMfgQ3UHt6yv283IXEVGrx+BDdAcRHOBMRNRmMPgQ3cGgm2d8jl0yobyy2snVEBHRL8HgQ3QH/t7t4O/dDharQMb5EmeXQ0REvwCDD1ED1N7dtY8vLCUiatUYfIgaoHaA84FzDD5ERK0Zgw9RA9QGn4PnSlBlsTq5GiIiaiwGH6IGuEfdAV7tXHGjyoITl0zOLoeIiBqJwYeoAeRyGQbfHOfD5/kQEbVeDD5EDVR7WzsHOBMRtV6NCj7JyckICgqCm5sbdDod9u3bd9v2JSUlMBgM0Gq1UKlUCA4ORkpKim357t278fjjj8PPzw8ymQxffPFFnT5eeOEFyGQyu2n06NF2bYqLizFx4kR4enrC29sbL774Iq5du9aYXSSq48cBzlchhHByNURE1BgOB581a9YgNjYW8fHxOHjwIMLCwjBq1CgUFhbW295sNmPEiBHIzc3FunXrkJWVhWXLlsHf39/Wpry8HGFhYUhOTr7ttkePHo38/Hzb9Nlnn9ktnzhxIo4fP47t27dj06ZN2L17N6ZNm+boLhLVq7+/F1QuchSXm3HmcrmzyyEiokZwcXSFhQsXYurUqZgyZQoAYMmSJdi8eTOWL1+OOXPm1Gm/fPlyFBcXY8+ePXB1dQUABAUF2bUZM2YMxowZc8dtq1Qq+Pr61rvs5MmT2Lp1K/bv349BgwYBAD744AM8+uij+Otf/wo/Pz9HdpOoDqWLHAMCvLH3bDEO5Bajl7qDs0siIiIHOXTGx2w2Iz09HZGRkT92IJcjMjISaWlp9a6zceNG6PV6GAwGaDQa9OvXDwkJCbBYLA4X+80330CtVqN3796IiorClStXbMvS0tLg7e1tCz0AEBkZCblcjr1799bbX2VlJUwmk91EdDu1l7v2cYAzEVGr5FDwKSoqgsVigUajsZuv0WhgNBrrXScnJwfr1q2DxWJBSkoK5s6di7/97W94++23HSp09OjR+Ne//oXU1FS8++672LVrF8aMGWMLUEajEWq12m4dFxcXdOrU6Za1JSYmwsvLyzYFBAQ4VBNJz2DbC0uvOrkSIiJqDIcvdTnKarVCrVZj6dKlUCgUCA8Px8WLF5GUlIT4+PgG9/PMM8/Yfu7fvz9CQ0PRs2dPfPPNNxg+fHijaouLi0NsbKzts8lkYvih27ov0BtyGZBXfB0FpgpoPN2cXRIRETnAoTM+Pj4+UCgUKCgosJtfUFBwy7E3Wq0WwcHBUCgUtnkhISEwGo0wm82NKLlGjx494OPjg+zsbACAr69vnQHW1dXVKC4uvmVtKpUKnp6edhPR7Xi4uSJEW/N7wuf5EBG1Pg4FH6VSifDwcKSmptrmWa1WpKamQq/X17vO0KFDkZ2dDav1x8f8nzp1ClqtFkqlspFlAxcuXMCVK1eg1WoBAHq9HiUlJUhPT7e12bFjB6xWK3Q6XaO3Q/RzteN89vN5PkRErY7Dt7PHxsZi2bJlWLlyJU6ePImoqCiUl5fb7vKaNGkS4uLibO2joqJQXFyM6OhonDp1Cps3b0ZCQgIMBoOtzbVr15CRkYGMjAwAwNmzZ5GRkYG8vDzb8tdffx0//PADcnNzkZqaiieffBK9evXCqFGjANScRRo9ejSmTp2Kffv24fvvv8eMGTPwzDPP8I4ualI/DnDmOB8iotbG4TE+EyZMwOXLlzFv3jwYjUYMGDAAW7dutQ14zsvLg1z+Y54KCAjAtm3bEBMTg9DQUPj7+yM6OhqzZ8+2tTlw4AAefvhh2+facTeTJ0/GihUroFAocOTIEaxcuRIlJSXw8/PDyJEj8dZbb0GlUtnWW7VqFWbMmIHhw4dDLpdj/PjxeP/99x0/KkS3UfvqikyjCaaKKni6uTq5IiIiaiiZ4CNobUwmE7y8vFBaWsrxPnRbw5J24tyV6/hkymA83Ft95xWIiKjZOPL9zXd1ETWC7fUVHOBMRNSqMPgQNULEzeCTvPOMkyshIiJHMPgQNYKuR03wcZHLUGiqcHI1RETUUAw+RI3QrXN73BfojWqrwNr0C84uh4iIGojBh6iRfhcRCABYvT8PVivvESAiag0YfIga6bFQP3i4ueB88Q18f6bI2eUQEVEDMPgQNVI7pQJPD/QHAHy2L8/J1RARUUMw+BD9Ar/T1Vzu+up4AS6XVTq5GiIiuhMGH6JfoI+vJwbaBjmfd3Y5RER0Bww+RL+QbZDzvvMc5ExEdJdj8CH6hR4L1cJD5YK84uvYc+aKs8shIqLbYPAh+oXclS4Yx0HOREStAoMPUROovdy17biRg5yJiO5iDD5ETaCvnyfCAmoGOf/vIJ/kTER0t2LwIWoiz0YEAABW7+OTnImI7lYMPkRN5PEwP3RQuSD3ynX8kMNBzkREdyMGH6ImUjPI2Q8A8CkHORMR3ZUYfIia0E8HOV+5xkHORER3GwYfoiZ0r58Xwrp6ocoisC6dg5yJiO42DD5ETaz2rM9n+/IgBAc5ExHdTRh8iJrY42F+aK9UIPfKdaRxkDMR0V2FwYeoibVXueBJ25Oc+eJSIqK7CYMPUTN4tnaQ8zEOciYiupsw+BA1g37+Xujv7wWzxYrPD150djlERHQTgw9RM6kd5PxOyklcLLnh5GqIiAhg8CFqNk8O8ENAp3YAgN8uSUNuUbmTKyIiIgYfombSXuWCNdP06O7THhdLbuA3H6Uhy1jm7LKIiCSNwYeoGfl5t8Oal4agj68HLpdVYsLSNBy9UOrssoiIJIvBh6iZqT3csHraEIR19ULJ9So8u+wH7M8tdnZZRESSxOBD1AK83ZX4z+91iOjeCWWV1Zj08T58e/qys8siIpIcBh+iFuLh5oqVUyIwLLgLblRZ8OKKA/jquNHZZRERSQqDD1ELaqdUYOmkcIy+1xdmixVRqw5iQwaf80NE1FIYfIhamMpFgcXPDsTTA/1hsQq8siYDQXM2O7ssIiJJYPAhcgIXhRx//U0YJuoCUfsC98QtJ2G18m3uRETNicGHyEnkchneHtcPsx7pBQD4aFcOolal47q52smVERG1XQw+RE4kk8kQO7I3/j4hDEqFHNuOF+C3H6XBWFrh7NKIiNqkRgWf5ORkBAUFwc3NDTqdDvv27btt+5KSEhgMBmi1WqhUKgQHByMlJcW2fPfu3Xj88cfh5+cHmUyGL774wm79qqoqzJ49G/3790f79u3h5+eHSZMm4dKlS3btgoKCIJPJ7KYFCxY0ZheJWtRTA7vi06k6dGqvxLGLJjyZ/B2OXeSDDomImprDwWfNmjWIjY1FfHw8Dh48iLCwMIwaNQqFhYX1tjebzRgxYgRyc3Oxbt06ZGVlYdmyZfD397e1KS8vR1hYGJKTk+vt4/r16zh48CDmzp2LgwcP4vPPP0dWVhaeeOKJOm3ffPNN5Ofn26aZM2c6uotETjEoqBO+mD4UvdQdUGCqxG+WpGEbb3cnImpSMiGEQ6MpdTodBg8ejMWLFwMArFYrAgICMHPmTMyZM6dO+yVLliApKQmZmZlwdXW9c0EyGdavX49x48bdtt3+/fsRERGBc+fOITCw5i3YQUFBeOWVV/DKK680aF8qKytRWVlp+2wymRAQEIDS0lJ4eno2qA+iplZ6owozPj2Ib08XQSYD5ozug2kP9oBMJnN2aUREdyWTyQQvL68GfX87dMbHbDYjPT0dkZGRP3YglyMyMhJpaWn1rrNx40bo9XoYDAZoNBr069cPCQkJsFgsjmy6jtLSUshkMnh7e9vNX7BgATp37oyBAwciKSkJ1dW3HiiamJgILy8v2xQQEPCLaiJqCl7tXPHJC4Px3JCaO74St2Ri9v+OwFxtdXZpREStnkPBp6ioCBaLBRqNxm6+RqOB0Vj/KfmcnBysW7cOFosFKSkpmDt3Lv72t7/h7bffbnTRFRUVmD17Nn73u9/ZJbtZs2Zh9erV2LlzJ1566SUkJCTgD3/4wy37iYuLQ2lpqW06f/58o2siakouCjneerIf4h/vC7kM+O+BCwj+8xZcuVZ555WJiOiWXJp7A1arFWq1GkuXLoVCoUB4eDguXryIpKQkxMfHO9xfVVUVfvvb30IIgQ8//NBuWWxsrO3n0NBQKJVKvPTSS0hMTIRKparTl0qlqnc+0d1AJpNhytDuCOrcHjM+PYhyswVPLP4eHz0fjn7+Xs4uj4ioVXLojI+Pjw8UCgUKCgrs5hcUFMDX17fedbRaLYKDg6FQKGzzQkJCYDQaYTabHSq2NvScO3cO27dvv+N1PJ1Oh+rqauTm5jq0HaK7ycN91PjCMBTdfdrjYskNjP9wD9YfuuDssoiIWiWHgo9SqUR4eDhSU1Nt86xWK1JTU6HX6+tdZ+jQocjOzobV+uP4hFOnTkGr1UKpVDZ427Wh5/Tp0/j666/RuXPnO66TkZEBuVwOtVrd4O0Q3Y3u0XjgC8NQPNy7CyqrrYhZcxhvbTqBagvH/RAROcLh29ljY2OxbNkyrFy5EidPnkRUVBTKy8sxZcoUAMCkSZMQFxdnax8VFYXi4mJER0fj1KlT2Lx5MxISEmAwGGxtrl27hoyMDGRkZAAAzp49i4yMDOTl5QGoCT2//vWvceDAAaxatQoWiwVGo9HurFFaWhoWLVqEw4cPIycnB6tWrUJMTAyee+45dOzYsdEHiOhu4dXOFf+cPBgzHq550vPH353FpOX7UFzu2JlTIiJJE43wwQcfiMDAQKFUKkVERIT44YcfbMuGDRsmJk+ebNd+z549QqfTCZVKJXr06CHeeecdUV1dbVu+c+dOAaDOVNvP2bNn610OQOzcuVMIIUR6errQ6XTCy8tLuLm5iZCQEJGQkCAqKioavF+lpaUCgCgtLW3MYSFqMVuOXhIhc7eIbrM3ifsTU8XRCyXOLomIyGkc+f52+Dk+bZkjzwEgcrYsYxmm/fsAzl25DjdXOd4dH4onB/jfeUUiojam2Z7jQ0R3j96+Htho+BWGBXdBRZUV0aszEDRnM5/3Q0R0Gww+RK2Yl7srlr8wGNMf6mmb95sle5B35boTqyIiunsx+BC1cgq5DH8Y3QcfPR8Or3auOHyhFGPf/xZfHr5055WJiCSGwYeojRh1ry9Soh/AoG4dUVZZjZmfHcKc/x3BDfMvez0MEVFbwuBD1Ib4e7fD6mlDMPORXpDJgNX7z+OJxd8h02hydmlERHcFBh+iNsZFIcerI3tj1Ys6dPFQ4XThNTy5+Hus2nsOvImTiKSOwYeojbq/lw+2RD+Ah24+7flP649h+qqDKL1R5ezSiIichsGHqA3z6aDC8smD8adHQ+Ail2HLMSPC5n+FlKP5PPtDRJLE4EPUxsnlMkx9sAf+F3U/unV2BwBMX3UQk5bvw9micidXR0TUshh8iCQiLMAb2155ELMe6QWlQo5vTxdh1N93429fZfHOLyKSDAYfIglxc1UgdmRvbIt5EA8Gd4HZYsUHO7Ix4u+78PWJAmeXR0TU7Bh8iCSou097rJwyGEueuw9+Xm64cPUGfv+vA/j9yv04X8ynPhNR28WXlP4EX1JKUnTdXI33U7Pxz29zUG0VULnIUVltxel3xsBVwb+NiOjux5eUElGDuStdMGdMH2x95QHoe3RG5c2XnD65+Hscu1jq5OqIiJoWgw8RAQB6qT3w6VQdFk0YgI7urjiRb8K45O/xt6+yUFnNwc9E1DYw+BCRjUwmw7iB/vgqZhjG9tei2irwwY5sPP7Bd8g4X+Ls8oiIfjEGHyKqo4uHCskT78OHE++DTwclThVcw9P/+B6JKSdRUcWzP0TUejH4ENEtjemvxfaYYRg3wA9WAXy0OwePvvctDuQWO7s0IqJG4V1dP8G7uohu7esTBfjTF0dRYKqETAYIAXw3+2F07eju7NKISOJ4VxcRNbnIvhp8FTMMvx3UFbV/Lg1L+gazPjvEu7+IqNVg8CGiBvNq54q//DoMn00dgl/18oHFKrDx8CU89sF3eO6fe7H71OUmffnp6YIyvPzvdATN2Yz80htN1i+1fkIIfHXciKA5m7E35wpfuksNxktdP8FLXUSOOXaxFMu+zcGmI/mwWGv+VxKi9cS0B7vjsVC/Rj8AsbjcjL9vP4VP9+XZ+u3cXonFz94Hfc/OTVY/tU5pZ67g3a2Zdnca9vBpj98ODsD4+7qii4fKecW1ERarwOnCMmTkleDwhRLkFl1HWs4V5C4Y6+zS6uXI9zeDz08w+BA1zvni61j+/Vms2X8e12++8NTPyw2XSiuQ+uow9OzSoUH9VFZb8K895/D+jtMoq6gGAESGaHCx5AZO5pugkMsQN6YPXvxVd8hksmbbH7o7Hb9Uir9szcKuU5cBAO1cFXjgHh98l11k+71zkcswPESNZwYH4sHgLlDI+XvSEPmlN3D4fAkOnS9BRl4Jjl4stR3TWkqFHP+YeB8i+2qcVOWtMfg0EoMP0S9Tct2M//xwDiv25KLomtk2v0eX9hjRV4ORfX0xMMAb8p99GQkhsO24EYlbMnHuSs27wkK0npj7WAju7+mDG2YL4j4/gi8yLgEAHg/zw7vj+8Nd6dJyO0dOc+5KOf721SlsPFzz7+8il+GZiADMeuQeqD3dcK2yGpsOX8Lq/eftzgJpvdzwm/Cu+M2gAAR04iD8+qQczcf0VQfrXdZeqUBoV28MCPTG6YIyfH2yEC5yGT743UCM6a9t4Upvj8GnkRh8iJpGRZUFGzIuYtORfKSduYJq64//m/HpoEJkiBoj+mowtJcPsguv4a1NJ7D3bM0t8l08VHh9ZG+MD+9q99e6EAIr9+Ti7c0nUW0V6K3xwEfPhyPIp32L7x81nrnaii3H8hG9OgOzHumFzh1U8Omggk8HJXw8an72dHOBTCbD5bJKfLDjND7dm2f7HXo8zA+vjgi+5b97ptGENfvPY/2hiyi5XmWbPyDAG4/298WYflqGoJs+25eHP64/CiEAhVyGYI0HBgR4Y2BATdjp2aWD7b/BaosVsf89jI2HL0Ehl2Hhb8Pw5AB/J+/Bjxh8GonBh6jpmSqq8E3WZWw/UYBvMgtRVlltW9bOVYEbNx+IqHKRY+oDPfDyQz3RQXXrMzn7c4sxfdVBXC6rhIebC957ZgAe6XP3nXone1fLzfh0Xx7+lZaLAlPlbdsqXeTwaa/EpdIK27wHg7vgD6N6o5+/V4O2V1FlwVcnCrBmfx6+z75it+xeP0+M6eeL0f206KVu2GXYtmbJrjNYsCUTAPCsLhB/HhtyxzOoFqvA7P8dwbr0C5DJgHfHh+K3gwJaotw7YvBpJAYfouZlrrZi79kr+Op4Ab4+WYD8m19sTw7wwx9G94G/d7sG9VNgqsD0VQeRfu4qAOCVyHsw65F76lxCI+fLLryG5d+fxecHL6CiquYFuF08VHg81A9miwVFZWYUXau8OZlx7SfBGADCArwxe3Rv3N/Tp9E1FJgqsO24EVuOGrH37BX85AQk7lF3wJh+vng0VIs+vm3///tCCPxlWxY+/OYMACDqoZ74w6jeDR4zZ7UK/OmLY/hsXx4A4J2n+mGirluz1dtQDD6NxOBD1HKEEDiZX4Z2SgW6N+Jylbnairc2ncC/fzhnm5f9zhi4NPJOMikTQuDDXWfwl61Z+O9LekR07/SL+/v2dBGWf38W32Rdts2/188TL/6q5o4/pUv9/04VVRZcLqvElfKaMWJhXb2adCD7lWuV2H6iAFuOGbHnTBGqLD9+BQ4L7oJXRwYjtKt3k23vbmK1CszdcAyr9taEltmj+yDqoZ4O9yOEwPwvT2DFnlwAQPzjfTFlaPemLNVhDD6NxOBD1PqsS7+AP60/ispqK/4wujemP9TL2SW1Kj//EgOA30UEYs6YPvBq5+pQX9UWK77IuITX1h62zZPJgBEhGrz4q+6I6N7prrobr/RGFVJP1oSgnZmFtnFEo+7VIHZEb/T29XByhbdmrrbi8IUS9NV6ov1tLg3XqrJY8erNMToyGfD2uF92pkYIgQVbM/HRrhwAQNyYPnhpmOMhqqkw+DQSgw9R67Qu/QJeW3sYSoUcKdG/Qi/13fuFdTexWAX+/MVRfLbvPICaMx61t4p38VDhjcfvxaP9fe8YVmoDzwc7TtvuymuvVOC3gwPwwv1B6Nb57h+Afu5KOd77+jTWZ1yEEDWB7YkwP7wSGdyoM5LNqcBUgZf+nY6M8yVQusjxQC8fjLrXF8ND1Ojcoe4zjCqqLJi+6iB2ZNbclbVwwgA8Eeb3i+sQQuDvX5/G+6mnAQCxI4Ixa/g9v7jfxmDwaSQGH6LWSQiBKSv245usyxgQ4I3/Rd3P57fcQbXFitfXHcH6QxchlwF/+XUYfh3eFXtzriBu/VHkXC4HAAzvo8ab4/rVO/6qvsDTqb0S0x7sgWd1gfB0c+yM0d3gdEEZ/v71KaQcNQKoudvp1/d1xczhvezeS3fdXI0zheU4VVCG04XXcLqgDKcKy3C++AYSn+6P30UENkt9h/Ku4qV/p6OwrBIKucz2gE8AkMuAQUGdMOpeX4y6V4OuHd1RVlGFF1cewL6zxVC5yLHkuXA83EfdpDUl78xG0rYs2+f0P0fWG8CaE4NPIzH4ELVe+aU3MHLhbpRVVuOPj/bBtAedd9r9ToQQSDtzBQVlFXiktwZe7i0bEMzVVkSvPoQtx4xwkcuw6JkBeCz0xzMAldUW/GPnGfzjm2xUWQTclQq8NrI3Jt8fBIVcdsvA89KDPfC8vlubeL7SsYulWLj9FHZkFgKoeXjfo/19YaqoxunCMly4egO3+/ac9mAPzBndp0kH3K9Lv4A/fn4UZosVwZoOWDZpECqqrNh23Ihtx404fslk1/5eP09UWwSyCsrgoXLBxy8M/sXjt27ln9/m4J2UkxAC8HZ3xR8fDcFvwru22KVNBp9GYvAhat3W7M/D7P8dhcpFjpToBxr8xOiWUlltwYaMS/j427PIKigDUHMb/5h+vngmIhC6FhgD89PLHkqFHIufHYiR9/rW2/Z0QRn+uP4o9ufW3D0X2tULTw30x4o9uW028Pxc+rli/O2rU9hz5kqdZZ3bK3GPpgOCNR64R90B92g8kHbmCt67eelnZF8NFj0z4Bcfl2qLFYlbMvHxd2cBACP6avD3CQPqPPbhfPF1fHWiANuOG3Egt9h291qn9kr86/8iGvwogMY6lHcVcZ8fRaax5ndb170T3nmqf4s8MoDBp5EYfIhaNyEEJi3fh29PFyG8W0f89yX9XXHJ62p5zROtV6adQ9G1mmfYtFcqoPVuh+zCa7Z2QZ3dMWFwIMaH+0Pt4dbkdVw3V2Pav9LxXXYRVC5yLJ00CMOCu9x2HatVYPX+80jcctL2GhGg7Qeen9uTXYTdp4vg7+2Ge24GnVtdztmQcRGvrz0Cs8WK/v5e+OfkQdB4Nu7fs+S6GTM+PYTvsosAALOG34NXht/50Q1XrlUi9WQhjl0qxQv3B6FHC/0RUGWx4pPvz+Lv20/jRpUFrgoZoh7qhekP9YSbq6LZtsvg00gMPkSt38WSGxi5cBfKzRbMfawvXvyV826zzbl8DR9/dxb/+8kzbLRebnjh/iA8ExEITzcXHLlQitX787Ax4xLKb74bSSGXYXgfNZ6JCMCD93SBAHCtohrXKqtRbq62/XytshrlldWQQYaundohqHN7+Hq61fulWFZRhf9bsR/7c6/CXanAx5MHO/TC10JTBd7afBKH8q7i+SHdJBN4GutAbjGm/TsdxeVmaL3c8PHkwejr59j3yqmCMkz91wGcu3Id7VwVWPjbsLvuVRG3cr74OuZtOIadNx9n0N2nPd4Z1w/392r885hup9mDT3JyMpKSkmA0GhEWFoYPPvgAERERt2xfUlKCP/3pT/j8889RXFyMbt26YdGiRXj00UcBALt370ZSUhLS09ORn5+P9evXY9y4cXZ9CCEQHx+PZcuWoaSkBEOHDsWHH36Ie+75cQR5cXExZs6ciS+//BJyuRzjx4/He++9hw4dGpZ0GXyI2oZVe8/hT+uPwc1Vjq3RD7b4ay0O5V1F8s5spGYW2saB9PP3xNQHeuDR/tp631pfXlmNzUfysXp/Hg7mlTR620oXOQI7uSOoszsCO7VHkI87Ajq6Y1HqaRw+XwIPNxesmBKB8G4dG70NaphzV8rxfyv248zlcrRXKvDBswMb/JTxr44bEbMmA+VmC7p2bIdlkwYhRNu6vpeEENhyzIg3Nh5HYVnNmc6n7/PHnx4NafLBz80afNasWYNJkyZhyZIl0Ol0WLRoEdauXYusrCyo1XVHipvNZgwdOhRqtRp//OMf4e/vj3PnzsHb2xthYWEAgC1btuD7779HeHg4nn766XqDz7vvvovExESsXLkS3bt3x9y5c3H06FGcOHECbm41pxDHjBmD/Px8fPTRR6iqqsKUKVMwePBgfPrppw3aNwYforZBCIGJ/9yLPWeuIKJ7J6yeOqRFnup85VolFmzJxNr0C7Z5kSFq/P6BHg6N3zlVUIY1+8/j84MXcPUn75tSucjh4eaC9ioXtFe6oIObCzqoXFBtFThffB3ni6/bvRft5zq6u+LfL+qafawH/aj0ehWiVqVjz5krkMuAeY/1xQs/e9hfyXUzTuSbkJlfhkyjCSfzy3D0YikAYEiPTvjHxHB0aq90RvlNwlRRhb9uy8K/fzhn+0Mgd8HYpt1GcwYfnU6HwYMHY/HixQAAq9WKgIAAzJw5E3PmzKnTfsmSJUhKSkJmZiZcXe9854JMJqsTfIQQ8PPzw6uvvorXXnsNAFBaWgqNRoMVK1bgmWeewcmTJ9G3b1/s378fgwYNAgBs3boVjz76KC5cuAA/v7rPLKisrERl5Y/vjDGZTAgICGDwIWoDzhdfx6hFu3HdbMH8J+7F5PuDmm1bFqvAp/vykLQ1E6ab42DG39cV0x/u+YsGWJurrbh8rRLtlQq0V7nUe6bop6otVlwqqcC54nLkXrmOc0XlOFd8HeeulKODygWJT4fe1Q/la6uqLFb8ef0xrDlQ87yk3w7qis4dVMjMrwk5RlNFvetN0nfD3Mf63vHfvbU4lHcVf1p/DPGP94WuR8MvszZEswUfs9kMd3d3rFu3zi6YTJ48GSUlJdiwYUOddR599FF06tQJ7u7u2LBhA7p06YJnn30Ws2fPhkJRd6BTfcEnJycHPXv2xKFDhzBgwADb/GHDhmHAgAF47733sHz5crz66qu4evWqbXl1dTXc3Nywdu1aPPXUU3W29cYbb2D+/Pl15jP4ELUN/0rLxbwNx+GuVGBr9IMI7Nz0b+U+fL4Eczccw5ELNX+h99V64q1x/XgpiewIIbB0dw4Sb74Y9OcCOrVDiK8n+mg90VfrgXv9vNrkW+StVtEsZ18dCT4OjUwrKiqCxWKBRmN/jVKj0SAzs/5/zJycHOzYsQMTJ05ESkoKsrOzMX36dFRVVSE+Pr5B2zUajbbt/Hy7tcuMRmOdS20uLi7o1KmTrc3PxcXFITY21va59owPEbUNz+m6YfORfOw9W4zZ/zuCVb/XNdn/dK+Wm/GXbVlYvT8PQgAebi54bWRvPDek211xJxndXWQyGV4aVnMGcPX+PKg93RDi64EQrSd6+3rAoxU+7LEx7oYXCTf7kHyr1Qq1Wo2lS5dCoVAgPDwcFy9eRFJSUoODT3NRqVRQqVr26ZJE1HLkchneHR+K0e/tRlrOFXy6Lw/PDfllb5K2WgX+e+A83t2aaRt/8/R9/ogbE4IuHvz/Cd1eZF8NIvs2bIAzNQ+Hgo+Pjw8UCgUKCgrs5hcUFMDXt/4HYGm1Wri6utpd1goJCYHRaITZbIZSeecBW7V9FxQUQKv98Va+goIC26UvX19fFBYW2q1XXV2N4uLiW9ZGRG1fkE97/GFUH7y56QT+/MUxnC4ow+8f6OHwZYQqixWbj+Rj6e4cnMiveUJub40H3hrXr9mehktETc+hEVNKpRLh4eFITU21zbNarUhNTYVer693naFDhyI7OxtWq9U279SpU9BqtQ0KPQDQvXt3+Pr62m3XZDJh7969tu3q9XqUlJQgPT3d1mbHjh2wWq3Q6XSO7CYRtTEv3B+EyJCav7JXpp3DQ3/9BjFrMpBpNN1hzZrn3yzbnYMH/7ITr6zJwIl8EzqoXPDnsSHYNOtXDD1ErYzDl7piY2MxefJkDBo0CBEREVi0aBHKy8sxZcoUAMCkSZPg7++PxMREAEBUVBQWL16M6OhozJw5E6dPn0ZCQgJmzZpl6/PatWvIzs62fT579iwyMjLQqVMnBAYGQiaT4ZVXXsHbb7+Ne+65x3Y7u5+fn20QdEhICEaPHo2pU6diyZIlqKqqwowZM/DMM8/Ue0cXEUmHXC7DsknhSDtzBR/uOoNvTxdh/aGLWH/oIh7po0bUQz0xOMg+wFwquYEVe3Lx2d48lFXW3Knl00GFF+7vhom6bujYim8vJpIyh4PPhAkTcPnyZcybNw9GoxEDBgzA1q1bbQOP8/LyIJf/eCIpICAA27ZtQ0xMDEJDQ+Hv74/o6GjMnj3b1ubAgQN4+OGHbZ9rBxxPnjwZK1asAAD84Q9/QHl5OaZNm4aSkhL86le/wtatW23P8AGAVatWYcaMGRg+fLjtAYbvv/++o7tIRG2QTCbD/b18cH8vHxy7WIoPd53BlqP52JFZiB2ZhQjv1hFRw3rC18sNH393Fl8evmR7Jk4vdQdMfaA7nhzg36yP3Sei5sdXVvwEH2BIJC1ni8qxdHcO/pd+AWaLtc7yIT06YdqDPfBQsPquuBuFiOrHd3U1EoMPkTQVmiqw/PtcrPrhHK5XWfBofy2mPtAdoV29nV0aETUAg08jMfgQSVtFlQVVFqtknqlC1FY02wMMiYjaMjdXBcfwELVxbeMFIEREREQNwOBDREREksHgQ0RERJLB4ENERESSweBDREREksHgQ0RERJLB4ENERESSweBDREREksHgQ0RERJLB4ENERESSweBDREREksHgQ0RERJLB4ENERESSweBDREREksHgQ0RERJLB4ENERESSweBDREREksHgQ0RERJLB4ENERESSweBDREREksHgQ0RERJLB4ENERESSweBDREREksHgQ0RERJLB4ENERESSweBDREREksHgQ0RERJLB4ENERESSweBDREREksHgQ0RERJLB4ENERESSweBDREREktGo4JOcnIygoCC4ublBp9Nh3759t21fUlICg8EArVYLlUqF4OBgpKSkNLjP3NxcyGSyeqe1a9fa2tW3fPXq1Y3ZRSIiImqDHA4+a9asQWxsLOLj43Hw4EGEhYVh1KhRKCwsrLe92WzGiBEjkJubi3Xr1iErKwvLli2Dv79/g/sMCAhAfn6+3TR//nx06NABY8aMsdveJ598Ytdu3Lhxju4iERERtVEyIYRwZAWdTofBgwdj8eLFAACr1YqAgADMnDkTc+bMqdN+yZIlSEpKQmZmJlxdXZukTwAYOHAg7rvvPnz88cc/7oxMhvXr1zc67JhMJnh5eaG0tBSenp6N6oOIiIhaliPf3w6d8TGbzUhPT0dkZOSPHcjliIyMRFpaWr3rbNy4EXq9HgaDARqNBv369UNCQgIsFkuj+0xPT0dGRgZefPHFOssMBgN8fHwQERGB5cuX43a5rrKyEiaTyW4iIiKitsvFkcZFRUWwWCzQaDR28zUaDTIzM+tdJycnBzt27MDEiRORkpKC7OxsTJ8+HVVVVYiPj29Unx9//DFCQkJw//33281/88038cgjj8Dd3R1fffUVpk+fjmvXrmHWrFn19pOYmIj58+c3dPeJiIiolXMo+DSG1WqFWq3G0qVLoVAoEB4ejosXLyIpKQnx8fEO93fjxg18+umnmDt3bp1lP503cOBAlJeXIykp6ZbBJy4uDrGxsbbPJpMJAQEBDtdERERErYNDl7p8fHygUChQUFBgN7+goAC+vr71rqPVahEcHAyFQmGbFxISAqPRCLPZ7HCf69atw/Xr1zFp0qQ71qvT6XDhwgVUVlbWu1ylUsHT09NuIiIiorbLoeCjVCoRHh6O1NRU2zyr1YrU1FTo9fp61xk6dCiys7NhtVpt806dOgWtVgulUulwnx9//DGeeOIJdOnS5Y71ZmRkoGPHjlCpVI7sJhEREbVRDl/qio2NxeTJkzFo0CBERERg0aJFKC8vx5QpUwAAkyZNgr+/PxITEwEAUVFRWLx4MaKjozFz5kycPn0aCQkJdpef7tRnrezsbOzevbvOM4AA4Msvv0RBQQGGDBkCNzc3bN++HQkJCXjttdcc3UUiIiJqoxwOPhMmTMDly5cxb948GI1GDBgwAFu3brUNTs7Ly4Nc/uOJpICAAGzbtg0xMTEIDQ2Fv78/oqOjMXv27Ab3WWv58uXo2rUrRo4cWacuV1dXJCcnIyYmBkII9OrVCwsXLsTUqVMd3UUiIiJqoxx+jk9bxuf4EBERtT7N9hwfIiIiotaMwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgkg8GHiIiIJIPBh4iIiCSDwYeIiIgko1HBJzk5GUFBQXBzc4NOp8O+fftu276kpAQGgwFarRYqlQrBwcFISUlxqM+HHnoIMpnMbnr55Zft2uTl5WHs2LFwd3eHWq3G66+/jurq6sbsIhEREbVBLo6usGbNGsTGxmLJkiXQ6XRYtGgRRo0ahaysLKjV6jrtzWYzRowYAbVajXXr1sHf3x/nzp2Dt7e3w31OnToVb775pu2zu7u77WeLxYKxY8fC19cXe/bsQX5+PiZNmgRXV1ckJCQ4uptERETUFgkHRURECIPBYPtssViEn5+fSExMrLf9hx9+KHr06CHMZvMv6nPYsGEiOjr6ln2kpKQIuVwujEaj3bY9PT1FZWVlQ3ZNlJaWCgCitLS0Qe2JiIjI+Rz5/nboUpfZbEZ6ejoiIyNt8+RyOSIjI5GWllbvOhs3boRer4fBYIBGo0G/fv2QkJAAi8XicJ+rVq2Cj48P+vXrh7i4OFy/ft22LC0tDf3794dGo7HNGzVqFEwmE44fP15vbZWVlTCZTHYTERERtV0OXeoqKiqCxWKxCxcAoNFokJmZWe86OTk52LFjByZOnIiUlBRkZ2dj+vTpqKqqQnx8fIP7fPbZZ9GtWzf4+fnhyJEjmD17NrKysvD5558DAIxGY7191C6rT2JiIubPn+/IISAiIqJWzOExPo6yWq1Qq9VYunQpFAoFwsPDcfHiRSQlJSE+Pr7B/UybNs32c//+/aHVajF8+HCcOXMGPXv2bFRtcXFxiI2NtX02mUwICAhoVF9ERER093Mo+Pj4+EChUKCgoMBufkFBAXx9fetdR6vVwtXVFQqFwjYvJCQERqMRZrO5UX0CgE6nAwBkZ2ejZ8+e8PX1rXMnWG2ft+pHpVJBpVLdchtERETUtjg0xkepVCI8PBypqam2eVarFampqdDr9fWuM3ToUGRnZ8NqtdrmnTp1ClqtFkqlslF9AkBGRgaAmmAFAHq9HkePHkVhYaGtzfbt2+Hp6Ym+ffs6sptERETUVjk6cnr16tVCpVKJFStWiBMnTohp06YJb29v291Uzz//vJgzZ46tfV5envDw8BAzZswQWVlZYtOmTUKtVou33367wX1mZ2eLN998Uxw4cECcPXtWbNiwQfTo0UM8+OCDtj6qq6tFv379xMiRI0VGRobYunWr6NKli4iLi2vwvvGuLiIiotbHke9vh8f4TJgwAZcvX8a8efNgNBoxYMAAbN261TaQOC8vD3L5jyeSAgICsG3bNsTExCA0NBT+/v6Ijo7G7NmzG9ynUqnE119/jUWLFqG8vBwBAQEYP348/vznP9v6UCgU2LRpE6KioqDX69G+fXtMnjzZ7rk/REREJG0yIYRwdhF3C5PJBC8vL5SWlsLT09PZ5RAREVEDOPL9zXd1ERERkWQw+BAREZFkMPgQERGRZDD4EBERkWQw+BAREZFkMPgQERGRZDD4EBERkWQw+BAREZFkMPgQERGRZDD4EBERkWQw+BAREZFkMPgQERGRZDD4EBERkWQw+BAREZFkMPgQERGRZDD4EBERkWQw+BAREZFkMPgQERGRZDD4EBERkWQw+BAREZFkMPgQERGRZDD4EBERkWQw+BAREZFkMPgQERGRZDD4EBERkWQw+BAREZFkMPgQERGRZDD4EBERkWQw+BAREZFkMPgQERGRZDD4EBERkWQw+BAREZFkMPgQERGRZDD4EBERkWQw+BAREZFkNCr4JCcnIygoCG5ubtDpdNi3b99t25eUlMBgMECr1UKlUiE4OBgpKSkN7rO4uBgzZ85E79690a5dOwQGBmLWrFkoLS2160Mmk9WZVq9e3ZhdJCIiojbIxdEV1qxZg9jYWCxZsgQ6nQ6LFi3CqFGjkJWVBbVaXae92WzGiBEjoFarsW7dOvj7++PcuXPw9vZucJ+XLl3CpUuX8Ne//hV9+/bFuXPn8PLLL+PSpUtYt26d3fY++eQTjB492vb5p9shIiIiaZMJIYQjK+h0OgwePBiLFy8GAFitVgQEBGDmzJmYM2dOnfZLlixBUlISMjMz4erq2iR9AsDatWvx3HPPoby8HC4uNflNJpNh/fr1GDdunCO7ZGMymeDl5YXS0lJ4eno2qg8iIiJqWY58fzt0qctsNiM9PR2RkZE/diCXIzIyEmlpafWus3HjRuj1ehgMBmg0GvTr1w8JCQmwWCyN7hOAbedqQ08tg8EAHx8fREREYPny5bhdrqusrITJZLKbiIiIqO1y6FJXUVERLBYLNBqN3XyNRoPMzMx618nJycGOHTswceJEpKSkIDs7G9OnT0dVVRXi4+Mb1WdRURHeeustTJs2zW7+m2++iUceeQTu7u746quvMH36dFy7dg2zZs2qt5/ExETMnz+/obtPRERErZzDY3wcZbVaoVarsXTpUigUCoSHh+PixYtISkpCfHy8w/2ZTCaMHTsWffv2xRtvvGG3bO7cubafBw4ciPLyciQlJd0y+MTFxSE2Ntau74CAAIdrIiIiotbBoUtdPj4+UCgUKCgosJtfUFAAX1/fetfRarUIDg6GQqGwzQsJCYHRaITZbHaoz7KyMowePRoeHh5Yv379LccM1dLpdLhw4QIqKyvrXa5SqeDp6Wk3ERERUdvlUPBRKpUIDw9HamqqbZ7VakVqair0en296wwdOhTZ2dmwWq22eadOnYJWq4VSqWxwnyaTCSNHjoRSqcTGjRvh5uZ2x3ozMjLQsWNHqFQqR3aTiIiI2iiHL3XFxsZi8uTJGDRoECIiIrBo0SKUl5djypQpAIBJkybB398fiYmJAICoqCgsXrwY0dHRmDlzJk6fPo2EhAS7y0936rM29Fy/fh3/+c9/7AYid+nSBQqFAl9++SUKCgowZMgQuLm5Yfv27UhISMBrr732iw8SERERtQ0OB58JEybg8uXLmDdvHoxGIwYMGICtW7faBifn5eVBLv/xRFJAQAC2bduGmJgYhIaGwt/fH9HR0Zg9e3aD+zx48CD27t0LAOjVq5ddPWfPnkVQUBBcXV2RnJyMmJgYCCHQq1cvLFy4EFOnTnX8qBAREVGb5PBzfNoyPseHiIio9Wm25/gQERERtWYMPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5EREQkGY0KPsnJyQgKCoKbmxt0Oh327dt32/YlJSUwGAzQarVQqVQIDg5GSkqKQ31WVFTAYDCgc+fO6NChA8aPH4+CggK7Nnl5eRg7dizc3d2hVqvx+uuvo7q6ujG7SERERG2Qw8FnzZo1iI2NRXx8PA4ePIiwsDCMGjUKhYWF9bY3m80YMWIEcnNzsW7dOmRlZWHZsmXw9/d3qM+YmBh8+eWXWLt2LXbt2oVLly7h6aefti23WCwYO3YszGYz9uzZg5UrV2LFihWYN2+eo7tIREREbZRMCCEcWUGn02Hw4MFYvHgxAMBqtSIgIAAzZ87EnDlz6rRfsmQJkpKSkJmZCVdX10b1WVpaii5duuDTTz/Fr3/9awBAZmYmQkJCkJaWhiFDhmDLli147LHHcOnSJWg0Gtu2Z8+ejcuXL0OpVNbZbmVlJSorK22fS0tLERgYiPPnz8PT09ORw0JEREROYjKZEBAQgJKSEnh5ed2+sXBAZWWlUCgUYv369XbzJ02aJJ544ol61xkzZoyYOHGimDp1qlCr1eLee+8V77zzjqiurm5wn6mpqQKAuHr1ql2bwMBAsXDhQiGEEHPnzhVhYWF2y3NycgQAcfDgwXpri4+PFwA4ceLEiRMnTm1gOn/+/B2zjAscUFRUBIvFYjujUkuj0SAzM7PedXJycrBjxw5MnDgRKSkpyM7OxvTp01FVVYX4+PgG9Wk0GqFUKuHt7V2njdFotLWpr4/aZfWJi4tDbGys7bPVakVxcTE6d+4MmUx2h6PhmNo0yrNJzYvHuWXwOLcMHueWwePccprrWAshUFZWBj8/vzu2dSj4NIbVaoVarcbSpUuhUCgQHh6OixcvIikpCfHx8c29+dtSqVRQqVR2834erpqap6cn/8NqATzOLYPHuWXwOLcMHueW0xzH+o6XuG5yaHCzj48PFApFnbupCgoK4OvrW+86Wq0WwcHBUCgUtnkhISEwGo0wm80N6tPX1xdmsxklJSW3bVNfH7XLiIiIiBwKPkqlEuHh4UhNTbXNs1qtSE1NhV6vr3edoUOHIjs7G1ar1Tbv1KlT0Gq1UCqVDeozPDwcrq6udm2ysrKQl5dna6PX63H06FG7O8G2b98OT09P9O3b15HdJCIiorbqjqOAfmb16tVCpVKJFStWiBMnTohp06YJb29vYTQahRBCPP/882LOnDm29nl5ecLDw0PMmDFDZGVliU2bNgm1Wi3efvvtBvcphBAvv/yyCAwMFDt27BAHDhwQer1e6PV62/Lq6mrRr18/MXLkSJGRkSG2bt0qunTpIuLi4hzdxWZRUVEh4uPjRUVFhbNLadN4nFsGj3PL4HFuGTzOLeduONYOBx8hhPjggw9EYGCgUCqVIiIiQvzwww+2ZcOGDROTJ0+2a79nzx6h0+mESqUSPXr0sLurqyF9CiHEjRs3xPTp00XHjh2Fu7u7eOqpp0R+fr5dm9zcXDFmzBjRrl074ePjI1599VVRVVXVmF0kIiKiNsjh5/gQERERtVZ8VxcRERFJBoMPERERSQaDDxEREUkGgw8RERFJBoNPC0hOTkZQUBDc3Nyg0+mwb98+Z5fU6u3evRuPP/44/Pz8IJPJ8MUXX9gtF0Jg3rx50Gq1aNeuHSIjI3H69GnnFNtKJSYmYvDgwfDw8IBarca4ceOQlZVl16aiogIGgwGdO3dGhw4dMH78+DoPEqU7+/DDDxEaGmp7mq1er8eWLVtsy3mcm96CBQsgk8nwyiuv2ObxODeNN954AzKZzG7q06ePbbmzjzODTzNbs2YNYmNjER8fj4MHDyIsLAyjRo2ye9AiOa68vBxhYWFITk6ud/lf/vIXvP/++1iyZAn27t2L9u3bY9SoUaioqGjhSluvXbt2wWAw4IcffsD27dtRVVWFkSNHory83NYmJiYGX375JdauXYtdu3bh0qVLePrpp51YdevUtWtXLFiwAOnp6Thw4AAeeeQRPPnkkzh+/DgAHuemtn//fnz00UcIDQ21m8/j3HTuvfde5Ofn26bvvvvOtszpx9nJt9O3eREREcJgMNg+WywW4efnJxITE51YVdsCQKxfv9722Wq1Cl9fX5GUlGSbV1JSIlQqlfjss8+cUGHbUFhYKACIXbt2CSFqjqmrq6tYu3atrc3JkycFAJGWluasMtuMjh07in/+8588zk2srKxM3HPPPWL79u1i2LBhIjo6WgjB3+emFB8fL8LCwupddjccZ57xaUZmsxnp6emIjIy0zZPL5YiMjERaWpoTK2vbzp49C6PRaHfcvby8oNPpeNx/gdLSUgBAp06dAADp6emoqqqyO859+vRBYGAgj/MvYLFYsHr1apSXl0Ov1/M4NzGDwYCxY8faHU+Av89N7fTp0/Dz80OPHj0wceJE5OXlAbg7jnOzv51dyoqKimCxWKDRaOzmazQaZGZmOqmqts9oNAJAvce9dhk5xmq14pVXXsHQoUPRr18/ADXHWalUwtvb264tj3PjHD16FHq9HhUVFejQoQPWr1+Pvn37IiMjg8e5iaxevRoHDx7E/v376yzj73PT0el0WLFiBXr37o38/HzMnz8fDzzwAI4dO3ZXHGcGHyK6I4PBgGPHjtldp6em1bt3b2RkZKC0tBTr1q3D5MmTsWvXLmeX1WacP38e0dHR2L59O9zc3JxdTps2ZswY28+hoaHQ6XTo1q0b/vvf/6Jdu3ZOrKwGL3U1Ix8fHygUijqj1QsKCuDr6+ukqtq+2mPL4940ZsyYgU2bNmHnzp3o2rWrbb6vry/MZjNKSkrs2vM4N45SqUSvXr0QHh6OxMREhIWF4b333uNxbiLp6ekoLCzEfffdBxcXF7i4uGDXrl14//334eLiAo1Gw+PcTLy9vREcHIzs7Oy74veZwacZKZVKhIeHIzU11TbParUiNTUVer3eiZW1bd27d4evr6/dcTeZTNi7dy+PuwOEEJgxYwbWr1+PHTt2oHv37nbLw8PD4erqanecs7KykJeXx+PcBKxWKyorK3mcm8jw4cNx9OhRZGRk2KZBgwZh4sSJtp95nJvHtWvXcObMGWi12rvj97lFhlBL2OrVq4VKpRIrVqwQJ06cENOmTRPe3t7CaDQ6u7RWraysTBw6dEgcOnRIABALFy4Uhw4dEufOnRNCCLFgwQLh7e0tNmzYII4cOSKefPJJ0b17d3Hjxg0nV956REVFCS8vL/HNN9+I/Px823T9+nVbm5dfflkEBgaKHTt2iAMHDgi9Xi/0er0Tq26d5syZI3bt2iXOnj0rjhw5IubMmSNkMpn46quvhBA8zs3lp3d1CcHj3FReffVV8c0334izZ8+K77//XkRGRgofHx9RWFgohHD+cWbwaQEffPCBCAwMFEqlUkRERIgffvjB2SW1ejt37hQA6kyTJ08WQtTc0j537lyh0WiESqUSw4cPF1lZWc4tupWp7/gCEJ988omtzY0bN8T06dNFx44dhbu7u3jqqadEfn6+84pupf7v//5PdOvWTSiVStGlSxcxfPhwW+gRgse5ufw8+PA4N40JEyYIrVYrlEql8Pf3FxMmTBDZ2dm25c4+zjIhhGiZc0tEREREzsUxPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5EREQkGQw+REREJBkMPkRERCQZDD5EREQkGf8PoHpfLuQp2DsAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(range(start, min(stop, start + len(plaquette))), plaquette, \",-\")\n",
    "plt.ylim(0.6, 0.62)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": ".venv",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.2"
  },
  "rise": {
   "scroll": true,
   "slideNumber": "c/t"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
